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during 4X scan: a) 2D image of the cross-section in XY plane, b1-b11) 2D images of the cross- 
Section alone the XZ, Plame oa te uesetasmueemnandae eee oeutesteasencrta sue gent) 165 
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Figure 5-19. (a) Photograph of Mantis Shrimp. (b) Schematic image of the transverse section of 
the dactyl club (highlighted in Figure 5-18(a)) where I, II, and II display the impact, periodic, and 
striated regions, respectively. (c) SEM micrograph of a fractured surface of the periodic region. 
(d-e) Representative prism creation: the multilayer section highlighted in multiple colors is built 
from the orientation of each layer with a pitch angle Y. .........sssosssssoeenssssssssceerssssssssecersssssssseeeeees 166 


Figure 5-20. Numerical framework representing the first step of multiscale modeling of the 
mechanical response of (a) a four-filament configuration with semi-rectangular shape in the 
filament axis, knowing the tensile strength of filament and interface (Tnfil and Tnint, exp) and 
fracture energy of filament and interface (G/cfil and GIcint, exp) obtained from experiment; (b) 
a four-filament configuration with rectangular-shape in the filament axis, knowing the tensile 
strength and fracture energy of filament (Tnfil and GIcfil) obtained from experiment and tensile 
strength and fracture energy of the new interface (Tnint, num and GIcint, num) calculated from 
numen CAN IMO C Ine a. cherie acaba cht ctr dedaiiad daieta aioe A 168 


Figure 5-21. Validation of cohesive properties of four-filament free of gap configuration with of 
those for four-filament configuration with diamond-shaped gaps in terms of averaged stress- 
displacement response. The insets represent the evolution of the cracking of interface at 
corresponding loading applied to the configuration. .......... cece cccccccceeseeeeseeceeeeeeeeeeeseeeeeeeeeeeaeeeeees 169 


Figure 5-22. Mechanical response of lamellar architecture using the three-point flexural test: a) 
flexural stress-displacement (6;-6) for unrestrained (undamaged) and restrained (damaged) printed 
prisms with characteristic angles of 0 = O and y = 0, as well as unrestrained printed, prisms with 
characteristic angles of 0 = 90 and y = 0. The insets represent the evolution of cracking at 
corresponding loading applied to the hcp element. 20.0... ccccccccesssseeeeeceeeeeeeaeesseeeeceeeeeaaaeenses 171 


Figure 5-23. Mechanical response of Bouligand architecture using the three-point flexural test: a) 
stress-displacement for unrestrained and restrained printed specimens; b) specific modulus of 
rupture (MOR); c) schematic visualization of twisting crack pattern evolution in unrestrained and 
restrained Bouligand architecture with pitch angle y = 3°. ......ccccsssssesessseeeeeeeeceeeeeeeeeeeeaaaaas 174 
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ABSTRACT 


The exceptional long-term performance of concrete 1s a primary reason that this material 
represents a significant portion of the construction industry. However, a portion of this 
construction material is prone to premature deterioration for multi-physical durability issues such 
as internal frost damage, restrained shrinkage damage, and aggregate susceptibility to fracture. 
Since each durability issue is associated with a unique damage mechanism, this study aims at 
investigating the underlying physical mechanisms individually by characterizing the mechanical 
and thermal properties development and indicating how each unique damage mechanism may 
compromise the properties development over the design life of the material. 

The first contribution of this work is on the characterization of thermal behavior of porous 
media (e.g., cement-based material) with a complex solid-fluid coupling subject to thermal cycling. 
By combining Young-Kelvin-Laplace equation with a computational heat transfer approach, we 
can calculate the contributions of (1) pore pressure development associated with solidification and 
melting of pore fluid, (11) pore size distribution, and (111) equilibrium phase diagram of multiple 
phase change materials, to the thermal response of porous mortar and concrete during 
freezing/thawing cycles. Our first finding indicates that the impact of pore size (and curvature) on 
freezing 1s relatively insignificant, while the effect of pore size is much more significant during 
melting. The fluid inside pores smaller than 5 nm (1.e., gel pores) has a relatively small contribution 
in the macroscopic freeze-thaw behavior of mortar specimens within the temperature range used 
in this study (.e., +24°C to -35°C). Our second finding shows that porous cementitious 
composites containing lightweight aggregates (LWAs) impregnated with an organic phase change 
material (PCM) as thermal energy storage (TES) agents have the significant capability of 
improving the freeze-thaw performance. We also find that the phase transitions associated with 
the freezing/melting of PCM occur gradually over a narrow temperature range (rather than an 
instantaneous event). The pore size effect of LWA on freezing and melting behavior of PCM is 
found to be relatively small. Through validation of simulation results with lab-scale experimental 
data, we then employ the model to investigate the effectiveness of PCMs with various transition 
temperatures on reducing the impact of freeze-thaw cycling within concrete pavements located in 


different regions of United States. 
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The second contribution of this work is on quantification of mechanical properties 
development of cementitious composites across multiple length scales, and two damage 
mechanisms associated with aggregate fracture and restrained shrinkage cracking that lead to 
compromising the long-term durability of the material. The former issue is addressed by 
combining finite element method-based numerical tools, computational homogenization 
techniques, and analytical methods, where we observe a competing fracture mechanism for early- 
age cracking at two length scales of mortar (meso-level) and concrete (macro-level). When the 
tensile strength of the cement paste is lower than the tensile strength of the aggregate phase, the 
crack propagates across the paste. When the tensile strength of the cement paste exceeds that of 
the aggregate, the cracks begin to deflect and propagate through the aggregates. As such, a critical 
degree of hydration (associated with a particular time) exists below which the cement paste phase 
is weaker than the aggregate phase at the onset of hydration. This has implications on the inference 
of kinetic based parameters from mechanical testing (e.g., activation energy). Next, we focus on 
digital fabrication of a cement paste structure with controlled architecture to allow for mitigating 
the intrinsic damage induced by inherent shrinkage behavior followed by extrinsic damage exerted 
by external loading. Our findings show that the interfaces between the printed filaments tend to 
behave as the first layer of protection by enabling the structure to accommodate the damage by 
deflecting the microcrack propagation into the stable configuration of interfaces fabricated 
between the filaments of first and second layers. This fracture behavior promotes the damage 
localization within the first layer (1.e., sacrificial layer), without sacrificing the overall strength of 
specimen by inhibiting the microcrack advancement into the neighboring layers, promoting a novel 
damage localization mechanism. This study is undertaken to characterize the shrinkage-induced 
internal damage in 7-day 3D-printed and cast specimens qualitatively using X-ray 
microtomography (wCT) technique in conjunction with multiple mechanical testing, and finite 
element numerical modeling. As the final step, the second layer of protection is introduced by 
offering an enhanced damage resistance property through employing bioinspired Bouligand 
architectures, promoting a damage delocalization mechanism throughout the specimen. This novel 
integration of damage localization-delocalization mechanisms allows the material to enhance its 
flaw tolerant properties and long-term durability characteristics, where the reduction in the 


modulus of rupture (MOR) of hardened cement paste (hcp) elements with restrained shrinkage 
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cracking has been significantly improved by ~ 25% when compared to their conventionally cast 


hcp counterparts. 
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1. INTRODUCTION 


1.1 Literature Review 


Cementitious composite materials are made of cement, water, and aggregate of various sizes. 
These materials are known for their exceptional performance over their designed service life. 
However, a considerable portion of the infrastructures made by these materials undergoes 
premature deterioration. According to the Department of Transportation, $231.4 billion per year 
is needed just to maintain the existing concrete pavement network, which has fallen into the 
backlog of roads in poor condition, in a good and operable condition. These durability-related 
issues are simultaneously governed by multiple complex damage mechanisms, depending on 
multiple primary factors, including intrinsic dynamic chemical reaction in concrete as well as 
extrinsic environmental conditions. It is noteworthy to mention that at a certain age of the specimen, 
a few specific mechanisms are responsible for being the primary reasons for premature 
deterioration, classifying this performance problem into two distinct categories of short-term and 
long-term durability issues. To this end, this study intends to provide a comprehensive multi- 
physical experimental-numerical framework to find the driving reasons that contribute to lowering 
the intrinsic resistance of cementitious materials over their service life period. The durability- 
related issues examined in this work are visually displayed in Figure 1-1. 

Cementitious composites exhibit a homogeneous continuum structure at the macro level; 
However, when these materials are scanned at lower length scales (across nano- to meso-level), 
they have a complex heterogeneous structure with phases of different physical properties. The 
primary goal of this study is to examine the principle components that play an important role in 


governing a specific durability issue. 
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Durability Issues in 
Cementitious Composites 





Shrinkage Aggrege Frost Damage 
Susceptibility to 
Fracture 


Figure 1-1. Diagram of the three common durability issues observed in cementitious systems. 


Probing the nano- and micro-levels, this study investigates the contribution of 
heterogeneities in governing the macroscopic physical behavior of cementitious systems across 
these length scales. Depending on the water/cement, type of aggregates, the addition of different 
admixtures, as well as environmental conditions (e.g., temperature and humidity), the rate of the 
chemical reaction between water and cement for making cement paste is determined. During the 
hydration process, the two phases of calcium silicate hydrate (C-S-H) and calcium hydroxide (CH) 
constitute the majority of products during cement hydration, where comprise about 50~60% and 
20-25% of the products in volume, respectively. The CH phase contributes to mainly control the 
pH of the solution at the level of 12.5. In contrast, the formation of C-S-H gel particles divided by 
thin layers of water (gel pores network) is mainly responsible for the gradual development of 
mechanical properties of cement-based materials. In simple words, the hardened cement paste is 
considered to include four main phases of unreacted cement clinker, CH, C-S-H gel (and its gel 
pores), and capillary pores. Figure 1-2 depicts these four main components. 

Due to a net decrease of solid product volume during hydration and water consumption, a 
portion of space of originally water-filled space between (cement) particles remains unfilled, 


known as capillary pore space. The disjoining pressure of solution in a capillary pore, arising from 
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dissimilar surface tensions of the solution and cement particle, leads to shrinking the C-S-H gel 
products, which causes a common type of durability-related issue in cementitious systems. In 
conjunction with the shrinkage issue, the internal frost damage is also a primary durability-related 
issue takes place from nanograins level to the mesoscale, where the phase transition of available 
solution entrapped within the pores of various sizes (e.g., gel pores, capillary pores, and air voids) 


leads to multiple degradation mechanisms. 
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Figure 1-2. Backscattered electron imaging of the hardened cement pastes at 1 day. Visualization 
of hydration productions of Calcium-Silicate-Hydrate (C-S-H), Calcium Hydroxide (CH), 
unhydrated residual cement particles, and pores of various sizes (adapted from 
Stutzman(Stutzman, 2000). 
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Cementitious composites, such as mortar and concrete, have complex structures where 
damage and failure mechanisms can occur across different length scales. The mechanical response 
of these materials is frequently predicted by considering the material structure as a homogeneous 
continuum at the macro-level. However, when these materials are viewed at the meso-level, they 
have heterogeneous structures in which concrete is composed of coarse aggregate inclusions 
surrounded by a mortar matrix formed by cement paste and fine aggregate inclusions (Folker H. 


Wittmann, 1983), shown in Figure 1-3. The key mesostructure morphological (1.e., size, shape, 
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volume fraction, and spatial distribution of aggregates), as well as mechanical properties (stiffness, 
strength, and fracture energy), play an important role in the observable mechanical macroscopic 
behavior. Prior meso-level numerical models were developed to analyze the heterogeneous and 
anisotropic behavior of concrete to obtain knowledge on the basic mechanisms underlying the 
observed macroscopic behavior (Du & Jin, 2014; Huang et al., 2015; Lopez et al., 2008a; Ren et 
al., 2015; B. Sun et al., 2015; X. Wang et al., 2015; X. F. Wang et al., 2015; A. Yin et al., 2015). 
While these mesomechanical modelings qualitatively investigated the overall mechanical response 
as well as fracture process in the mesostructure constituents including matrix, aggregate inclusions, 
and the interface transition zone between them (ITZ), little is known about the interplay between 
the evolution of hydrating cement paste phase at early ages, the morphological and the mechanical 
properties of aggregate phase, and how this combination affects the likelihood of early-age 


cracking in cementitious composites. 
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Figure 1-3. Microstructure of the interfacial region between cement paste and fine aggregate 
known as interfacial transition zone (ITZ), adapted from (Simonova et al., 2017). The scale bar 
indicates a length of 5 um. 
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1.1.1 Internal Frost Damage 


Cementitious composites are prone to cracking internally when they experience thermal 
cycling that enables freezing and thawing actions. The damage mechanisms associated with freeze- 
thaw actions include (1) the hydraulic pressure (1.e., negative tensile stress) attributed to the 
crystallization of in-pore solution which was firstly proposed by Powers (Powers, 1949), (2) the 
osmotic pressure induced by cryo-suction effect where initially frozen water in larger pores tend 
to absorb excess water from adjacent smaller unfrozen pores (Coussy & Monteiro, 2008; Powers 
et al., 1953), and (3) the ice crystallization pressure exerted by a thin layer of water between ice 
crystals formed in the pore and the pore wall (L. Liu et al., 2011). All these mechanisms are 
factually governed by the effect of the role of curvature (also known as pore size effect, or Gibbs- 
Thomson effect). The general idea is that as the temperature goes below freezing point, pore 
solution solidification begins inside the pores of larger size. Sun and Scherer (Z. Sun & Scherer, 
2010a) have performed extensive research to understand better the role of air voids (which have 
the largest size) on the physical behavior of mortar specimens experiencing internal frost damage. 
They measured the volume change of mortar under one freezing cycle, the amount of ice produced 
at specific stages of the freezing process, and the physical properties of frozen specimens for strain 
and stress calculations. In frozen saturated specimens without air-entrainment, the experiment 
observations showed that the hydraulic pressure is the primary mechanism responsible for damage. 

In contrast, ice crystallization pressure exerted to the mesopore walls was responsible for 
damage in air-entrained saturated specimens through transporting the considerable amount of 
solution from the mesopores to the air voids. Liu et al. (L. Liu et al., 2011) have carried out 
thermodynamic- and fracture-based simulations to explain the role of crystallization pressure as 
the primary source of frost damage in saturated cement paste at a micro-level. Similar to results 
reported by Sun and Scherer, these results confirmed the direct relation between volumetric 
dilation and an increasing number of microcracks with decreasing the temperature below the 
freezing point. Coussy and Monteiro (Coussy & Monteiro, 2008) provided a novel poromechanical 
approach to investigate the role of pore size distribution and suggested a theoretical formula for a 
spacing factor between the air voids connected through pore network channel for frost damage 
mitigation. This model, however, was incapable of capturing the micro-crack development since 
a homogenization scheme was employed to convert the pore pressure associated with freezing to 


the effective stress tensor. In order to accommodate this issue, Zhou et al. (T. Zhou et al., 2019) 
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developed a novel framework based on the distribution of pore solution and stress field at the 
nanograin level to indicate that the largest freezing-induced pressure is associated with the 
capillary forces concentrated at the interface between large capillary pores and gel pores. 
Consequently, a homogenized length scale for the representative volume element (RVE) was 
obtained for upscaling the physical behavior from nano- to macro-level. 

Based on various internal and external conditions, including the degree of saturation of 
specimen, mixture proportioning, rate of thermal loading, and size of the specimen, there are plenty 
of competing mechanisms to account for simulation of the physical behavior of cementitious 
composites with a complex liquid-solid coupling. Nevertheless, there is little known about at which 
temperature the ice nucleation begins to take place, the size range of freezing pores, and the volume 
of freezing pores at the freezing onset. Li et al. (W. Li et al., 2012) used a passive acoustic emission 
technique to monitor the progress of internal frost damage of mortar specimens with 13% air 
content and 96% degree of saturation. After performing two consecutive thermal cycles, a 
significant number of acoustic events were observed once the temperature dropped below -8°C, 
highlighting the high density of microcrack formation due to frost damage. Farnam et al. (Yaghoob 
Farnam, Bentz, Sakulich, et al., 2014) have provided sufficient evidence on both thermal and 
mechanical behavior of partially and fully saturated mortar specimens with air entrainment. 
However, there is still little known about the influence of the role of curvature and associated 


capillary pressure on the macroscopic physical behavior of porous cementitious systems. 


1.1.2 Frost Resistance Improvement Using Sustainable Organic Materials 


After analyzing the role of curvature on the pore solution solidification at the microscopic 
scale and assess the validity of resulting macroscopic thermal behavior of frozen mortar specimens 
with experimental data, we leverage from these earlier findings to investigate the effectiveness of 
sustainable green materials incorporated into cementitious composites on enhancing their internal 
frost resistance. A research study by Riayzi et al. (Riyazi et al., 2017) emphasized the use of 
superabsorbent polymers (SAPs) for synthetic air void manufacturing that leads to frost resistance 
enhancement of concrete. Organic (green) phase change materials (PCMs) can be used as 
substitutes to inorganic admixtures for thermal energy storage (TES) purposes, since these 
materials are capable of recurrently storing and releasing latent heat energy upon solidification and 


melting, respectively. Owing to their significant thermal mass, PCMs with a high transition 
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temperature (about 40 ~ 50 °C) were incorporated into concrete pavements to mitigate the thermal 
cracking by diminishing the interior temperature rise in early-age concrete (Dale P. Bentz & Turpin, 
2007; Fernandes et al., 2014). Bentz and Turpin also developed a computer simulation to highlight 
the effectiveness of PCM use on frost damage mitigation in bridge structures. Sakulich and Bentz 
(Aaron R. Sakulich & Bentz, 2012) provided evidence on the feasibility of low-temperature PCM 
addition for service-life extension of bridge decks in various climates located throughout the 
United States. 

To (1) incorporate and retain the PCMs in the cementitious composites, (11) indulge any 
potential chemical reactions between the fresh cement paste and PCM, and (111) maximize the 
thermal efficiency of cementitious composites by optimizing the distribution of PCMs, two 
practical methods of micro-encapsulation and lightweight aggregate (LWA) impregnation were 
used. Two recent research studies by Yeon and Kim (Yeon & Kim, 2018) and Urgessa et al. 
(Urgessa et al., 2019) have proven the effectiveness of low-temperature PCMs (with a 
characteristic freezing point of 4.5 °C) microencapsulated with melamine-formaldehyde resin in 
improving the frost resistance of mortar and concrete specimens. Farnam et al. (Yaghoob Farnam 
et al., 2017) and Li et al. (W. Li et al., 2019) have demonstrated the significant heat storage capacity 
of PCM-LWA composite on melting surface ice and snow and delaying internal ice formation, 
respectively. In the conventional method of concrete casting, it is unlikely to control the 
distribution of particulates during the mixing process as these inclusions were randomly distributed. 
While there is no carefully designed spatial pattern of PCM capsulate distribution, the literature 
falls short of addressing how to govern the temporal distribution of latent heat release and storage 
rate during freezing and thawing cycles. Leveraging from our findings on the role of curvature 
effect on the thermal behavior of freezing saturated mortar specimens, prudently robust design of 
PCM-LWA composite allows us to maintain the internal temperature of the specimen above the 
freezing point of the solution by (1) selecting PCMs with two favorable physical properties of 
latent heat energy and thermal diffusivity, in conjunction with (2) selecting a porous LWA 


particulate with a pore size range similar to that for cementitious composite. 


25 


1.1.3 Aggregate-matrix composite effect on the mesoscale fracture of early-age 
cementitious composites 


The prediction of early-age physical properties development of cementitious composites is 
of importance in various infrastructures, such as bridges and pavements. To this end, an empirical 
approach of maturity method is used to establish a direct relationship between strength 
development of material, cement hydration progress, and environmental conditions, including 
ambient moisture and curing temperature (Carino, 1991). The prediction accuracy of this method, 
however, is subjected to debate due to a lack of deep understanding of the relationship between 
hydration kinetics, microstructure development, maturity, and physical properties development 
(Jieying Zhang et al., 2008). As a result, Zhang et al. (Jieying Zhang et al., 2008) proposed a 
mathematical methodology to back-calculate the maturity control variable (1.e., activation energy) 
by lessening the error between empirical and experimental data. While this research highlighted 
that there is no unique value of activation energy for every physical property of concrete, we 
question the validity of this research finding as the activation energy merely depends on the 
mixture proportioning of the cementitious material and is independent of the type of external 
loading. 

On the other hand, this method analyzes the structure mortar and concrete as a 
homogeneous continuum, where the role of mesostructure phases of aggregate inclusions and ITZ 
on resulting physical properties is neglected. Based on the concept of structure-property, there is 
a strong effect of existing heterogeneous characteristics of mesostructure (both morphological and 
physical) on the resulting homogenized physical response of cementitious composites (1.e., size 
effect) (Elias et al., 2015). Outstanding works by Carpinteri (Carpinteri, 1986), and Bazant (Bazant, 
1997) have suggested a bottom-up approach for fracture simulation, explicitly considering a 
heterogeneous mesostructure. 

Previous investigations demonstrated the role that aggregates played in the early-age 
flexural strength development, highlighting a common short-term durability issue of concrete 
pavements (Barde et al., 2005). The aggregate volume content and the type of aggregate were 
varied to illustrate that the aggregates play an important role in strength development. Figure 1-4(a) 
shows a linear relationship between the flexural strength (f-) development and non-evaporable 
water (Wnevap) (1.e., assumed to be directly proportional to the quantity of chemically bound water) 


for cement paste specimen (0% aggregate). In comparison, the mortar (paste and 55% fine 
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aggregate by volume) and concrete (paste, 34% of fine aggregate and 38% of coarse aggregate by 
volume) show a bilinear response where aggregates are present. Also, two concrete specimens 
containing only 30% of a single size coarse aggregate were tested. Figure 1-4(b) illustrates how 
the failure strength is influenced by the aggregate type (e.g., limestone 1s weaker than granite) after 


the knee point, which implied that the strength of aggregate may influence the maturity prediction. 
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Figure 1-4. Experimental measurement of (a) effect of aggregate content on early age flexural 
strength development of cementitious composites, (b) effect of aggregate strength on flexural 
strength development of concrete, and photos of fractured surfaces of concrete specimens (c) at 
t = 12 h, and (d) t= 72 h (Barde et al., 2005). 
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The knee point is determined through a two-step process: first, the data points are clustered 
into two groups that will define the two linear rate of strength development. This is attained by 
maximizing each coefficient of determination (R^). Then the knee point is selected as the 
interception between the two linear regression functions. Barde et al. (Barde et al., 2005) also 
observed that the ratio of the number of fractured aggregates to the total number of aggregates in 
the fractured surface, defined as aggregate fracture density (AFD), significantly increased after the 
knee point of the bilinear behavior, illustrated by dashed lines in Figure 1-4(a) and (b). These 
findings imply that aggregate fracture strongly influences the knee point. Figure 1-4(c) and (d) 
show an example of two fractured surfaces of concrete beams subjected to four-point flexural test 
at a very early age of t = 12 h (corresponding to Wnevay = 14) and a later age of t=72h 
(corresponding tO Wnevap = 16.05), respectively. Figure 1-4(c) shows that the majority of coarse 
ageregates remained undamaged “and pulled out of the matrix,” and the crack propagates around 
the aggregates (before the knee-point). After the knee point, however, the aggregate fracture 
density suggests that the cracks have more tendency to propagate across the aggregates, shown in 


Figure 1-4(d). 


1.1.4 Mitigation of Shrinkage Cracking Through Digital Fabrication Approach 


Shrinkage of cementitious composites is an inevitable major durability issue which is 
caused due to the migration of water from the interior structure (1.e., higher water concentration) 
to the outer surface of the specimen (1.e., lower water concentration), and the net decrease of solid 
products volume during the hydration process. A variety of factors, including mixture 
proportioning, ambient conditions, curing age, the specimen size, and the restraint degree, are the 
key factors that determine the development of residual tensile stresses (S. Shah et al., 1998). When 
there is significant shrinkage, the cementitious composite undergoes a high density of 
microcracking development, which results in reducing the mechanical capacity of the specimen. 
Therefore, shrinkage-induced damage occurs at very early ages, and thus may significantly 
compromise the durability and service-life span of the cementitious composites (Holt & Leivo, 
2004). 

Five major shrinkage mechanisms are responsible for the volume reduction of cementitious 
composites: 1) chemical shrinkage occurs in a fresh mixture before hardening stage when the rate 


of water evaporation exceeds that of the bleeding water moving towards the surface of specimen, 
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2) drying shrinkage that refers to the capillary tension induced by in-pore condensed solution at 
the pore wall during moisture exchange with environment, causing some dimensional change of 
solid structure (Dale P Bentz et al., 1998), 3) autogenous shrinkage represents the volume 
contraction of cementitious composters in isothermal and sealed condition due to self-desiccation 
(generating empty pores within cement paste structure during hydration), 4) thermal shrinkage 
which corresponds to the volume contraction of massive structures due to significant rate of 
hydration and associated temperature development in hardening cementitious composites, and 5) 
carbonation shrinkage happens as a reaction between CH dissolution and ambient carbon dioxide 
on the surface of specimen and has no contribution in internal microcrack development. While 
thermal and carbonation shrinkage mechanisms can occur in concrete structures, these two 
shrinkage mechanisms have an insignificant impact on internal microcrack development in thin 
structures (such as pavements and bridges). Conversely, three mechanisms of plastic, autogenous, 
and drying shrinkage are the primary cause of early-age shrinkage in fresh and hardening cement- 
based materials (D. P. Bentz & Jensen, 2004; Cohen et al., 1990; Ghourchian et al., 2019; Lura et 
al., 2003), which all highlight the effect of water movement in volumetric contraction in 
conventional cast specimens. 

The current literature has surveyed multiple solutions for reducing the effect of shrinkage 
cracking on microcrack development and resulting macroscopic behavior of cementitious 
composites manufactured by conventional casting techniques. Weiss and Shah (W. J. Weiss & 
Shah, 2002) have shown the effectiveness of shrinkage reducing admixtures (SRAs) in reducing 
the early-age drying shrinkage by lowering the surface tension of the solution. Ghourchian et al. 
(Ghourchian et al., 2018a) provided a comprehensive understanding of the effect of chemical 
admixture incorporation on plastic shrinkage mitigation by modifying the material property rather 
than altering environment conditions (1.e., evaporation rate). Among various SRA (stabilizer), 
viscosity modifying admixture (VMA) (stabilizer), and C-S-H seeds, SRA was found to be the 
most efficient admixture by adjourning the capillary pressure rise. At the same time, no substantial 
impact on mechanical properties development during early stages of hydration was observed (no 
retarding effect). Radlinska et al. (Radlinska et al., 2008) clarified that SRA is a desirable 
admixture for shrinkage mitigation in a sealed system. In contrast, lightweight aggregate (LWA) 
incorporation lowers the shrinkage by providing sufficient amount of water distributed throughout 


the specimen volumetrically. 
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While the academic community has extensively explored several shrinkage mitigation 
strategies for the traditionally cast specimens, little research has been conducted on the capacity 
of digitally fabricated (1.e., 3D printed) cementitious composites on compensating the cracking 
associated with shrinkage. Buswell et al. (Buswell et al., 2018) have anticipated that the 3D printed 
concrete specimens are more likely to be susceptible to premature cracking due to significant 
autogenous shrinkage and greater exposed surface area when exposed to a dry environment. 
Marchon et al. (Marchon et al., 2018) have introduced the plastic shrinkage mechanism as the 
primary source of cracking development, wherein the absence of framework, excessive water loss 
may retard the hydration and resulting physical properties gain. Similarly, De Schutter et al. (De 
Schutter et al., 2018) observed the severe shrinkage cracking in a showcase of 3D printed cement 
mortar at an early age. To accommodate for creation of early-age cracks associated with drying 
and plastic shrinkage, Kazemian et al. (Kazemian et al., 2017) have used Polypropylene fiber as a 
reinforcement mechanism, given that the incorporation of this fiber of high tensile strength can 
significantly improve the tensile cracking resistance of a printed structure. Slavcheva (Slavcheva, 
2019) has provided quantitative evidence on increased deformability of fabricated cement paste 
structure using a 3D printing technique compared to that measured in the conventional casting 
technique. As a result, the magnitude of internal tensile stress in 3D printed hardened cement paste 
(hcp) elements may substantially increase as the drying shrinkage mechanisms is of the primary 
source of development of residual stresses under cyclic wetting-drying conditions. However, prior 
studies have failed to quantitatively investigate the adverse impact of shrinkage cracking on 
physical properties of 3D printed hcp elements and compare the results with those for traditionally 


cast hcp elements. 


1.2 Motivations for this study 


While cementitious composites interact with their surrounding environment, they frequently 
experience extensive physical changes that are likely to have a substantial adverse influence on 
their intrinsic properties. Thus, the durability of cementitious composites and their constituent 
phases at multiple length scales have received extensive consideration; however, several short- 
and long-term durability issues are of important concern that yet remained to be fully addressed. 
The level of complexity of these issues becomes more immense when the onset of each 


deterioration may begin at an unsightly lower length scale, which makes it more difficult to 
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understand the underlying physical damage mechanisms. This study will, first, focus on the role 
of curvature of pores of nano- to millimeter size on macroscopic freeze-thaw behavior of a 


cementitious composite with a complex fluid-solid coupling. 


Motivation I: To understand the effect of the role of curvature on macroscopic freeze-thaw 
behavior of air-entrained cementitious composites under thermal loading. 


Prior research has thoroughly investigated various thermodynamic- and physical-based 
mechanisms that are responsible for phase transition of the in-pore solution, solution transport 
within the porous network, and resulting volumetric change at the local and global scale. These 
studies proffered the damage onset associated with the freezing of solution within air voids in fully 
saturated cementitious composites, following the Gibbs-Thomson effect (Z. Sun & Scherer, 
2010a). However, little research has been conducted to identify a synergistic effect between the 
phase transition of the in-pore solution and the undercooling effect before ice formation onset. 
Owing to the undercooling phenomenon, in-pore ice crystals nucleate at a temperature lower than 
the characteristic freezing point of unconfined pore solution. The onset of freezing action, thus, 
may occur within multiple pore-length scales, including capillary pores (>5 nm of radius size) and 
air voids (>5 um). We, then, hypothesize that an instantaneous freezing action takes place in a 
broad range of pore sizes which promotes the contribution of capillary pores into frost damage 
onset, as observed using passive acoustic emission experiments performed by Farnam et al. 


(Yaghoob Farnam, Bentz, Sakulich, et al., 2014). 


Motivation 2: To understand the effect of sustainable phase change materials as thermal energy 
storage (TES) agent on macroscopic freeze-thaw behavior of cementitious composites. 


Deicing salts are often used to melt the ice or snow from the surface of concrete sidewalks 
and pavements to reduce the risk of injuries for the traveling public (W. Li et al., 2012). The water- 
salt solution that is produced during melting can infiltrate into the concrete and cause frost damage 
as the temperature drops below the freezing point of the solution in the pores (Yaghoob Farnam, 
Dick, et al., 2015; Yaghoob Farnam, Wiese, et al., 2015; Litvan, 1976; Powers, 1945; G. W. 
Scherer & Valenza, 2005; Z. Sun & Scherer, 2010c). Phase change materials (PCMs) are green- 
sustainable substances with a high enthalpy of fusion, 4H, that can be used to increase the thermal 


inertia of pavements(Cabeza et al., 2011). Hence, PCMs can be incorporated as alternatives to 
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deicing salts in some applications to reduce the number of freeze-thaw cycles experienced by 
pavements (Dale P. Bentz & Turpin, 2007; Yaghoob Farnam et al., 2017). Reducing freeze-thaw 
cycles is important as this may improve the mechanical response of pavement, thereby reducing 
the cracking. The presence of cracks can exacerbate reinforcing steel corrosion by providing 
aggressive media easy access to the steel and amplifying the amount of damage caused by 
successive freeze-thaw cycles (Kan et al., 2010). 

The main characteristic of a PCM 1s to release a considerable amount of thermal energy 
during the exothermic solidification reaction that occurs when the concrete pavement temperature 
drops to the freezing temperature of PCM. If the solidification reaction occurs at a slightly higher 
temperature than the freezing point of pore solution in concrete, it can delay the freezing of the 
pore solution (Yaghoob Farnam, Krafcik, et al., 2015; Aaron R. Sakulich & Bentz, 2012). The 
released heat maintains the surrounding pore solution at the freezing temperature of PCM, where 
the length of time that this temperature can be maintained is dictated by 4H; of PCM(Baetens et 
al., 2010). Many PCMs have high 44; and associated narrow temperature range for releasing and 
absorbing the thermal energy from exothermic solidification or endothermic melting reactions, 
respectively(He et al., 2004). This range corresponds to the phase transition temperature of the 
PCM. While the effectiveness of using PCM to delay or prevent freeze-thaw cycling is recognized 
as a promising technology, previous research has failed to investigate the contribution of the rate 
of 4H; and associated temperature range. 

Given the fact that LWA represents a porous aggregate material containing an internal 
porous network with a significant capacity of liquid absorption (Balapour & Farnam, 2018; Castro, 
Keiser, et al., 2011), we leverage our findings of a robust analysis on the pore size effect on the 
macroscopic thermal behavior of freezing and melting cementitious composite, to include two 
important thermodynamic-based effects of pore size distribution inside LWAs in conjunction with 


AH,;-temperature relationship obtained from equilibrium phase diagram of PCM. 


Motivation 3: To understand the intertwined effect between maturity method and composite nature 
of cementitious composites on their early-age fracture behavior. 


In the transportation environment, one primary factor of short-term durability issue is the 
premature failure of concrete pavements before reaching their designed physical properties (e.g., 


strength). The maturity method is used as an empirical concept that establishes an early-age 
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relationship between maturity index (1.e., degree of hydration which represents the fraction of 
cement that has chemically reacted with water) and strength of concrete (Carino, 1991). However, 
the contradicting field observations of premature deterioration of concrete pavements, particularly 
observed in high-performance concrete (HPC) pavements, highlights that the maturity method falls 
short in the prediction of accurate early-age properties. Therefore, more research is needed to 
understand better the competing fracture mechanism, which governs the resulting physical 
properties and time of cracking. 

One of the major drawbacks of the maturity method is modeling of cementitious 
composites as homogeneous isotropic materials, and merely considers the effects of ambient 
conditions on hydration process and properties development of cement paste. However, when 
these materials are viewed at one length scale lower than macro-level (1.e., meso-level), they 
exhibit heterogeneous structures where concrete is comprised of aggregates of larger sizes 
bounded by a mortar matrix which is developed by mixture of aggregates of smaller sizes and 
cement paste matrix (Folker H. Wittmann, 1983). While previous mesomechanical numerical tools 
examined the effect of mesostructure morphological properties as well mechanical properties on 
fracture behavior in the heterogeneous mortar and concrete, prior research has disregarded the 
strong interaction between the progress of cement paste hydration in conjunction with the 
morphological and mechanical properties of aggregates, and the effect of this combination on the 


probability of premature cracking in cementitious composites. 


Motivation 4: To gain insight into the synergistic effect between heterogeneous (macro-)structure, 
provided by the interface, and carefully designed architecture inspired by nature on shrinkage- 
induced damage in additively manufactured cement-based materials. 


As discussed in motivation 3, precise estimation of early-age physical properties of HPC is 
imperative in many structural applications such as concrete pavements, slabs, and bridges (W. J. 
Weiss et al., 2000). Besides the incorrect application of maturity method, a primary short-term 
durability issue of these thin (1.e., low thickness-to-width and low thickness-to-length ratios) 
structures with low w/c are premature cracking associated with multiple shrinkage mechanisms. 
Three water-related shrinkage mechanisms of plastic, autogenous, and drying are most likely to be 
responsible for the development of in-pore tensile stresses, which leads to detrimental early-age 


micro- and macro-cracking. Given the fact that prior research has extensively explored the 
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underlying physical mechanisms responsible for shrinkage cracking (Cohen et al., 1990; 
Ghourchian et al., 2018b; Lura et al., 2003; William Jason Weiss, 1999), and thus provided 
mitigation strategies for alleviating the corresponding detrimental effect on short-term durability 
of traditionally concrete cast specimens (D. P. Bentz & Jensen, 2004; Ghourchian et al., 2018a; 
Radlinska et al., 2008), the rise of novel processing technique of digital fabrication (also known 
as additive manufacturing and 3D printing) allows for controlling the architecture and associated 
anisotropy behavior of extrudable cementitious composites (Biernacki et al., 2017; Gosselin et al., 
2016; Lim et al., 2012; Perrot & Amziane, 2019; Sanjayan & Nematollahi, 2019). As a result, the 
fabrication of 3D printed specimens creates a unique opportunity by utilizing the creation of 
interface regions between filament layers and thus the creation of an extra level of anisotropy, 
which leads to arresting the shrinkage associated microcrack development and preventing their 
adverse effect on the mechanical performance of specimen. In simple words, while printed 
filament behaves identical to their cast counterparts under shrinkage loading, the creation of 
interface regions as “cold spots” for shrinkage-induced microcrack growth results in localizing the 
progressive damage. The claim of this novel proposal is in contrary with the current literature, as 
it is considered that fast rate of shrinkage deformation in filament layers, combined with 
mechanical properties dissimilarity of neighboring printed filaments can result in significant 
deformation and causing damage in the 3D printed specimen, ultimately causing a quicker 


premature failure compared to cast specimens. 


1.3 Objectives of this study 


This study aims at establishing a fundamental understanding of complex physical behavior of 
cementitious composites at multiple length scales, and extending our knowledge to (1) develop 
novel sustainable infrastructure materials that exhibit paradigm-shifting properties, (2) provide 
robust guidelines on enhancing their service life through combining the state-of-the-art 
computational and experimental findings to distinguish the primary physical mechanisms that lead 
to a specific durability issue. As a result of a comprehensive analysis of the physical behavior of 
damaged specimens, appropriate solutions to reducing the risk of premature cracking are proposed. 
The sub-objectives are found below: 

1) To identify factors that affect the internal frost damage in porous cementitious composites 


with inherent fluid-solid coupling. 
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2) 


3) 


4) 


To harness the incorporation of phase change material for alleviating the frost damage in 
porous cementitious composites. 

To understand the primary competing fracture mechanism in cementitious composites at 
multiple length scales, and its effect on the prediction of tensile mechanical properties by 
maturity (empirical) method. 

To learn about an intertwined relationship between the presence of heterogeneous 
interfaces fabricated by additive manufacturing, and novel bio-inspired architecture, how 
these two structural features promote a novel and complex localization-delocalization 


damage mechanism. 


1.4 Organization of Thesis 


This thesis contains five chapters. Its layout is addressed as follows: 


Chapter |: Introduction 

Chapter 2: Computational heat transfer modeling of wet cement mortar: the role of 
curvature effect 

Chapter 3: Enhancing the freeze-thaw performance of cementitious composites through the 
incorporation of phase change material (PCM) 

Chapter 4: A two-step multiscale model to predict early-age strength development of 
cementitious composites considering competing fracture mechanisms 

Chapter 5: Toward fabrication of enhancing damage-tolerant architecture cement-based 
material 


Chapter 6: Conclusion 


In the first chapter, the literature review, motivations, and objectives of this work are 


described thoroughly. Besides, the structure and organization are deliberated in the first chapter. 

In the second chapter, the interplay between the role of curvature and associated capillary pressure, 
and the undercooling phenomenon on the macroscopic thermal response of air-entrained mortar is 
investigated. This chapter simulates the heat flow and associated phase transformation events 


during thermal cycling in partially and fully saturated mortar specimens. This chapter was used to 
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develop a paper published in 2017 titled “Numerical simulation of the freeze-thaw behavior of 
mortar containing deicing salt solution” (Hadi S. Esmaeeli et al., 2017); and the contributions of 
Prof. Amir Yaghoob Farnam, Mr. Dale Bentz, Prof. Pablo Zavattieri, and Prof. Jason Weiss are 
greatly acknowledged. 

In the third chapter, the thermal performance of organic paraffin oil as a substitute for 
deicing salt on reducing the number of thermal cycles and thus improving the freeze-thaw 
performance of a small-scale mortar specimen and a large-scale concrete specimen is analyzed. 
Next, the effectiveness of multiple phase change materials with various characteristic transition 
temperatures was examined. This chapter was used to develop a paper published in 2018 titled 
“Numerical analysis of the freeze-thaw performance of cementitious composites that contain phase 
change material (PCM)” (Hadi S. Esmaeeli et al., 2018). The contributions of Prof. Amir Yaghoob 
Farnam, Prof. John Haddock, Prof. Pablo Zavattier1, and Prof. Jason Weiss are greatly 
acknowledged. 

In chapter four, a mesoscopic mechanistic approach was used to investigate the composite 
effect of cementitious composites on the development of early-age mechanical properties. We 
developed a two-step homogenization numerical framework to enable the upscaling of the physical 
properties of heterogenous mortar. This chapter was used to develop a paper published in 2019 
titled “A two-step multiscale model to predict early-age strength development of cementitious 
composites considering competing fracture mechanisms” (Hadi S Esmaeeli et al., 2019); and the 
contributions of Dr. Mehdi Shishehbor, Prof. Jason Weiss, and Prof. Pablo Zavattieri are greatly 
acknowledged. 

Chapter five describes a novel damage mechanism in the cementitious system by 
modification of hcp elements at the architectural scale. We manifested the role of weaker 
interfacial regions on localizing the cracking induced by shrinkage loading followed by 
delocalizing the cracking exerted by mechanical loading, which yields to enhancing the 
mechanical properties of the hcp element. This chapter will be used to develop a paper, and the 
contributions of Dr. Mohammadreza Moin, Prof. Jan Olek, and Prof. Jason Weiss, and Prof. Pablo 
Zavattierl are greatly acknowledged. 

In chapter six, the results obtained throughout the research are summarized, and the summary of 


the conclusion of the research is provided. 
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2. COMPUTATIONAL HEAT TRANSFER MODELING OF WET 
CEMENT MORTAR: THE ROLE OF CURVATURE EFFECT 


This chapter contains work that was originally published in Springer Nature as “Hadi S. 
Esmaeeli, Yaghoob Farnam, Dale P. Bentz, Pablo D. Zavattieri1, W. Jason Weiss. Numerical 
simulation of the freeze-thaw behavior of mortar containing deicing salt solution. In Journal of 
Materials and Structures. Volume 50., pp. 1-20. Springer, Cham, 2017.” The original article has 
been used with permission as stated below. 

Reprinted by permission from Springer Nature: Springer, Cham, In: Journal of Materials and 
Structures. “Numerical simulation of the freeze-thaw behavior of mortar containing deicing salt 
solution.” Hadi S. Esmaeeli, Yaghoob Farnam, Dale P. Bentz, Pablo D. Zavattieri, W. Jason Weiss, 
COPYRIGHT (2017). 


2.1 Introduction 


Prediction of phase transformation within the pores requires an understanding of heat flow 
within a mortar during a freezing/thawing cycle (Han et al., 2006; Radjy, 1968). For this discussion, 
the term “latent heat” is used to denote the amount of energy released or absorbed during a phase 
transformation (formation or melting of ice or eutectic solid). The latent heat produced by the 
phase transformation of the pore solution can be used to quantify the amount of pore solution in 
concrete that freezes (Powers, 1945). Two main phenomena affect the freezing of pore solution in 
a mortar/concrete: (1) its pore size distribution and (2) undercooling. The former influences its 
freezing. Concrete pores are typically categorized into three main classes: 1) gel pores with a radius 
smaller than 5 nm that is associated with the formation of cement (binder) hydration products, 2) 
capillary pores that are the remnants of the original water-filled space between (cement) particles 
and commonly range from 5 nm to 5 um in radius, and 3) pores (voids) associated with entrained 
or entrapped air that range from 5 um to 10 mm (Castro, Bentz, et al., 2011; Kumar & 
Bhattacharjee, 2003; Whiting & Nagi, 1998; Young, 1988). The size of the pores in the concrete 
can influence the temperature at which freezing occurs. This is described using the Gibbs- 
Thomson equation (Z. Sun & Scherer, 2010a). A large fraction of water associated with pore sizes 
greater than 5 nm (1.e., capillary pores or pores associated with entrained or entrapped air) is 
susceptible to freezing at a temperature above -10 °C (Cai & Liu, 1998; Yaghoob Farnam, Bentz, 
Sakulich, et al., 2014; W. Li et al., 2012). According to the Gibbs-Thomson equation, the water 


absorbed in the gel pores will not begin to freeze until the temperature of the specimen drops to 


a1 


about -13 °C (Cai & Liu, 1998; W. Li et al., 2011; Z. Sun & Scherer, 2010b). It is also worth 
mentioning that the solution inside concrete pores (1.e., pore solution) contains different ionic 
species (such as Na*, K*, Ca**, and OH’) (Andersson et al., 1989) that depress its freezing 
temperature (Mehta, P. K. and Monteiro, 2006). The absorption of salt solution into the pores can 
further depress the freezing temperature of this pore solution, due to the presence of additional 
ions such as Cl (Beddoe & Setzer, 1988; Pigeon & Pleau, 2010; G. Scherer, 1999). 

Undercooling also influences freezing in concrete. While it is expected that a solution freezes 
at its characteristic melting point temperature, Tm, freezing usually occurs at a temperature (1.e., Tp 
lower than Tm. This reduction in freezing temperature is known as undercooling (Askeland & 
Pradeep, 2003; Debenedetti, 2003; Wilding, 1992) and is primarily since solidification (in most 
cases) requires the presence of or formation of nuclei that can trigger the freezing action. Once the 
heterogeneous nuclei are present in the liquid phase, ice crystals begin to nucleate/grow. 
Consequently, the latent heat of fusion is released into the undercooled liquid, increasing the 
temperature of the liquid toward Tm. The growth of ice continues until the temperature of the liquid 
reaches Tm (Askeland & Pradeep, 2003). Afterward, the temperature of the liquid remains at Tm 
until the entire liquid solidifies; this is known as thermal arrest (Askeland & Pradeep, 2003). After 
thermal arrest, the amount of ice increases gradually as the further temperature decreases. Melting, 
however, occurs gradually in the pores as the temperature of each pore reaches its Tm value 
(Yaghoob Farnam, Bentz, Sakulich, et al., 2014; Qian et al., 2014; Z. Sun & Scherer, 2010b). The 
amount of ice transformed to solution increases gradually as each set of larger pores, in turn, 
reaches the associated Tm (according to the Gibbs-Thomson equation) (Z. Sun & Scherer, 2010b, 
2010c). 

These two phenomena (1.e., pore size and undercooling) affect the freezing behavior of pore 
solution simultaneously in the mortar specimen and it may be essential to consider both in the 
simulation during freezing. In melting, however, only pore size influences the thawing behavior. 
Numerically, it is feasible to develop a theoretical model based on the heat transfer formulation to 
predict and simulate phase transformation and heat transfer in materials (L. C. Thomas, n.d.; Velraj 
et al., 1999). One-dimensional finite difference (D. Bentz, 2000; Dale P. Bentz & Turpin, 2007; 
Costa et al., 1998; Lecomte & Mayer, 1985; Velraj et al., 1999), the two-dimensional finite 
difference (Dale P. Bentz & Turpin, 2007; Costa et al., 1998; Ismail & Abugderah, 2000; Velraj 
et al., 1999; Zivkovic & Fujii, 2001), control volume (Fukai et al., 2003; Hamada et al., 2003; 
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Ismail & Da Silva, 2003; Voller & Swaminathan, 1993), and finite element (Gong & Mujumdar, 
1997; Rubinsky & Cravahlo, 1981; Yoo & Rubinsky, 2007) based methods all have been used to 
simulate such heat transfer problems. In this paper, a one-dimensional finite-difference model 1s 
used. In particular, this numerical method approximates the complex solid-liquid interactions in 
the porous mortar using a fixed grid method (Hibbert et al., 1988; Shamsundar & Sparrow, 1975). 
The computational model is applied to estimate the thermal behavior of mortar containing in-pore 
water under freeze-thaw cycles. The formation of ice is quantified by calculating the volume 
fraction of ice that is produced (Callister & Rethwisch, 2009; Smith & Hashemi, 2006). As the ice 


grows, the liquid to solid phase transformation releases latent heat, AH p , that increases the 


temperature of the material locally and slows down the ice growth. An empirical approach is used 
to account for the sudden latent heat release produced by undercooling. This model is also used to 
describe, analyze and interpret the experimental data obtained from low-temperature longitudinal 
guarded comparative calorimeter (LGCC) tests (Yaghoob Farnam, Bentz, Hampton, et al., 2014; 
Yaghoob Farnam, Bentz, Sakulich, et al., 2014). 


2.2 Numerical Simulation of Heat Transfer Coupled with Phase Transformation 


The main objective of this study is to predict the thermal response of a mortar (considered 
at a macroscopic scale) that is experiencing phase transformations during a reduction and 
subsequent increase in specimen temperature (1.e., a freezing and thawing cycle). The goal of the 
simulations is to quantify the fraction of pore solution that can freeze in an undercooled mortar 
specimen. The temperature of the specimen can be tracked by solving the heat (energy balance) 
equation and considering the frozen fraction of the pore solution. The governing equation for the 
heat transfer within a mortar specimen can be described using the energy balance equation 2-1 
(Incropera et al., 2007). 


O OT (x,t) OT (x,t) 
“Ík (T) -q = p, (T\)-CP (T) _ 
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where T (x,t) is the temperature at location x(mm) and time t(sec), km (T) is the thermal 


conductivity of the mortar specimen [W/m K)| at temperature T, Pm (T) is the density of the 


mortar specimen (ks / m`) at temperature T, CË is its specific heat capacity [J i (kg-K)| at 
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temperature T, dye, 18 the rate of generated or consumed heat from any phase change of the pore 


solution E / (m sec) and Qjo., 1s the rate of heat dissipation (to the environment) in the 


experiment E / (m ‘sec )| ; 


In Equation 2-2a, the incorporation of a released/absorbed latent heat term, qoen associated 


with the freezing/melting of the pore solution within a mortar specimen is described. A heat sink 
term jos, 18 also included, as shown in equation 2-2b, to calculate the rate of heat dissipation to 
the environment (even though insulation is present). This heat term is considered as a fraction of 


the rate of generated latent heat to simulate the significant heat exchange between the mortar 


specimen and its surroundings in the lateral directions. 


OWT )- AT 
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where AH ; is the latent heat of fusion (kJ / kg) ; Psoln 18 the density of pore solution (kg / m’), 


Vp 1s the total volume fraction of pores within the mortar specimen (0 to 1), vp (T) is the volume 
fraction of the pore solution that can freeze at temperature T (0 to 1), é (T) is the frozen volume 


fraction of freezable pore solution with salt at temperature T (0 to 1), hoss is the heat dissipation 


coefficient (< 1), and AH f = AH ¢ -(1— hoss) is the apparent latent heat measured considering 


heat dissipation during phase transformation in the system (< AA). 


2.3 Frozen fraction of pore solution without salt, v-(T) 


The mortar specimens contain a pore structure with a broad range of sizes. The pore size can 
alter the freezing temperature of water (Beddoe & Setzer, 1988; Kaufmann, 2000; Mehta, P. K. 


and Monteiro, 2006; G. Scherer, 1999). To determine the pore size distribution in the mortar 
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specimen and thus to calculate v p (T) , a desorption isotherm was obtained for the mortar using a 


dynamic vapor sorption analyzer (TA Q5000). The vacuum saturation method is used to fully 
saturate the mortar specimen (1.e., Ds = 100 %). Therefore, all of the pores, including air voids, are 
filled with water, to investigate the role of curvature of the pores on the thermal behavior of the 
mortar. For melting, the pore size distribution obtained from an absorption isotherm is used (Litvan, 


1972). Figure 2-1 provides the desorption-absorption isotherm for the mortar specimen, and it 


correlates the degree of saturation (D, ) to the relative humidity (RH ), which is the amount of 


water vapor present in the specimen expressed as a percentage of the amount needed for saturation 
at the same temperature (Villani et al., 2015). A characteristic hysteresis is observed in the 
absorption/desorption isotherm in Figure 2-1, at least partially due to the presence of “ink-bottle” 


pores. 


To calculate vp (T) , two approaches were evaluated in this study: (1) a model with explicit 


consideration of a continuous pore size distribution, and (2) a phenomenological model with 


consideration of only a discrete pore size distribution. The first approach considers the effect of all 


pore sizes on the freezing process and vp (T) varies continuously as the temperature changes. In 


the second approach, the effect of a discrete pore size distribution on ice formation inside the 
mortar specimen is simplified. A phenomenological model is adapted to simulate the freezing 
process of water inside the mortar specimen (it considers only two classes of pores- large pores 
that include all pores except gel pores (the capillary and air-entrained/entrapped pores) and small 


pores (known as gel pores)). 


In the phenomenological model, vp (T) is considered to be a constant value based on three 


main classes of pores: 1) gel pores, 2) capillary pores, and 3) water-filled pores associated with 
entrained or entrapped air. To investigate the accuracy of these two approaches, the LGCC test 
conducted by Farnam et al. (Yaghoob Farnam, Bentz, Hampton, et al., 2014; Yaghoob Farnam, 
Bentz, Sakulich, et al., 2014) was simulated using these two models with consideration of 
continuous and discrete pore sizes, respectively. Also, the thermal behavior of the mortar 
specimens saturated with water was compared with experimental results obtained in a temperature 


range between 24 °C and -35 °C. 
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Figure 2-1. Desorption-absorption isotherms of mortar specimen. 


2.3.1 A model with consideration of a continuous pore size distribution 


Equation 2-3 describes the Gibbs-Thomson equation that relates the freezing temperature 


of a liquid inside a porous material to the pore radius. 


“FEL ii =c Kn -1;") = 
V 
r L 


where ycr, is the crystal/liquid interfacial energy (J / m°), r` is the radius of the pore for 
homogeneous nucleation (m), Sz and Sc are the molar entropies of the liquid and crystalline 
phases |J/(mol-K)|, Vr is the molar volume of the liquid (m°/mot , Ip is the melting 


temperature (K), and I’, is the freezing temperature as a function of pore radius (K)(Brun et al., 
1977; Z. Sun & Scherer, 2010b). Therefore, the temperature at which ice begins to form can be 
predicted as a function of the critical pore radius r by solving Equation 2-3 for Ty (Z. Sun & 


Scherer, 2010a, 2010b, 2010c). 


Figure 2-2 displays the relationship between the size of the pore and the temperature that 


is needed to freeze water inside the pore i: (°°). At a temperature above Ir (>); no phase 
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e e e e e X 
transformation occurs inside the pores with a radius less than r . Once the temperature reaches 


the associated freezing temperature, ice begins to form inside the pores with radii r . Afterward, 
ice propagates into the smaller pores, but only as the temperature drops further. This process is 


reversed during melting. 


Solution 
D d 


T=T.r) T>T(r) 





Figure 2-2. The effect of pore size on the freezing temperature of the water using the Gibbs- 
Thomson equation, including a schematic of ice formation in a porous material (inset). 


The Kelvin- Young-Laplace equation can be used to correlate the pore radius to the relative 


humidity (RH ) in a water-filled pore as described in Equation 2-4 (Henkensiefken et al., 2009; 


Radlinska et al., 2008). 


(aaah ; 
In( RH) ) \ RT 


In this study, the Kelvin- Young-Laplace equation (Equation 2-4) was used alongside the 





Gibbs-Thomson equation (Equation 2-3) to obtain the relationship between vp (T) and pore size 


in the mortar specimen based on its measured desorption isotherm (Figure 2-1). At a temperature 
of -35 °C, solution absorbed into the mortar pores with sizes greater than 1.47 nm is susceptible to 
freezing, as shown in Figure 2-2. Figure 2-2 also displays the process of ice formation in a porous 


material, as the ice forms inside the larger pores initially. Ice invades into the smaller pores 
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progressively as the temperature drops. Figure 2-3 shows that 72 % of the solution absorbed in the 
pores by volume can freeze between 0 °C and -35 °C. For the case of melting, the formed ice in 


the pores is similarly considered to melt gradually according to the Gibbs-Thomson equation 
(Equation 2-3). In this work, it is assumed that Il=v,,+(D,/100) and 
(D /100) = v, (7 < r*)+ Vr (r > r“) where Vair and v, are the volume fraction of air and non- 
freezable pore solution in the total pores (0 to 1), respectively. To investigate the role of pore sizes, 
all of the pores with various sizes are assumed to be filled with water in this section, i.e., V,;- =Q. 


It should be mentioned that the LGCC test was conducted within a temperature range between 


24 °C and -35 °C. 
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Figure 2-3. Volume fraction of pore solution that can freeze as a function of the critical nucleus 


(pore) size. 


2.3.2 A phenomenological model with consideration of a discrete pore size distribution 


Although the continuous pore size distribution can be estimated to determine the volume 
fraction of the freezable pore solution, this measurement is generally not available for 
field/commercial concretes. In this paper, a phenomenological model is developed to use in 
practice, when knowledge of the continuous pore size distribution is not available. Since the 
freezing and thawing responses of most cementitious systems are dominated by the category of 


relatively large pores (1.e., capillary, air-entrained, and air-entrapped pores), a discrete pore size of 
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5 nm (a pre-defined critical pore size as the division between gel and capillary pores) can be 
utilized as a criterion to differentiate the freezable and non-freezable pore solution (W. Li et al., 


2012; Z. Yang et al., 2006; Young, 1988). For a mortar specimen saturated with deionized water, 


this critical pore radius corresponds to a relative humidity (RH ) equal to 81 %. At this point, all 


of the gel pores are filled by a solution (Radlinska et al., 2008). 
The corresponding freezing temperature of water in pores with a size of 5 nm is about -13 


°C according to the Gibbs-Thomson equation as displayed in Figure 2-2. Figure 2-3 displays the 


relation between derived volume fraction of freezable pore solution, vp from associated relative 


. . . . . X . . . 
humidity measured in experiments, and pore radius, r , using Equation 2-4. The corresponding 
volume fraction of freezable pore solution, vr, in pores with size greater than 5 nm is measured to 


be 60 % concerning the total volume fraction of pores displayed in this figure. Consequently, the 
volume of freezable pore solution in the phenomenological model is underestimated by about 16 % 
concerning the model with direct consideration of pore size distribution, assuming pure water to 
be the solution in the pores. This implies that 60 % of the total solution by volume, corresponding 
to the solution that is absorbed into large pores (1.e., capillary pores, air-entrained pores, and air- 
entrapped pores) begins to transform to ice instantaneously. The gradual process of ice formation 
in the remainder of the freezable pores (i.e., smaller pores, containing a lower volume fraction) 
will be neglected. Therefore, the radius of curvature (pore size) would have a relatively small 
impact on the macroscopic freezing response of the air-entrained mortar specimen, and the 
approach of a discrete pore size distribution will be implemented in the numerical model to 
investigate the thermal behavior of mortar specimens containing NaCl solutions. 

To develop a reliable numerical simulation at the macro-scale, it 1s essential to define the material 
properties accurately. For detailed information on the utilization of homogenization techniques for 
calculation of the thermal properties of the mortar specimen, the readers are suggested to read my 


master’s thesis (Hadi Shagerdi Esmaeeli, 2015). 


2.4 Configuration of numerical simulation and boundary conditions 


Following the work by Farnam et al. (Yaghoob Farnam, Bentz, Sakulich, et al., 2014), the 
LGCC test was simulated to quantify heat flow and predict the temperature profiles of the mortar 


specimens. Two types of experimental data are used: 1) fully saturated mortar specimens (1.e., 100 % 
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degree of saturation) with solutions containing 0 %, 5 %, 10 %, and 23.3 % NaCl (by mass), and 
2) specimens saturated partially with water (1.e., no NaCl involved in the solution) at degrees of 
saturation equal to 75 %, 85 %, 95 %, and 100 %. The procedures used for the preparation of fully 
saturated and partially saturated mortar specimens were addressed in previous experimental works 
(Yaghoob Farnam, Bentz, Sakulich, et al., 2014; Yaghoob Farnam, Esmaeeli, et al., 2015; 
Yaghoob Farnam, Todak, et al., 2015). 

The experimental conditions of one-dimensional heat transfer were provided in the LGCC 
experiment by using a heat sink at the bottom, longitudinal insulation on the sides, and foam as 
thermal insulation around the system to minimize the heat dissipation from the experimental 
apparatus (Figure 2-4a). However, a difference between the measured released heat in an LGCC 
experiment and the associated enthalpy of fusion of phase change materials (1.e., Methyl Laurate 
and Paraffin Oil), likely due to experimental imperfections (thermal bridges, heat leaks, etc.), was 
observed. Therefore, Aioss, a heat loss coefficient, is employed to simulate the energy dissipation in 
the experimental system, which is estimated as a 40 % to 60 % heat loss (Yaghoob Farnam, Krafcik, 
et al., 2015). It is worth mentioning that the advection of heat to simulate the water transport 
occurring during the freeze/thaw cycle is neglected in this numerical investigation. 

Two references (meter) bars made of Pyroceram code 9606! with known thermal properties were 
used to measure the heat flow passing through the mortar specimen in the experiment (see Figure 


2-4a. 


' Certain commercial products are identified in this paper to specify the materials used and procedures employed. In no case does 
such identification imply endorsement or recommendation by the National Institute of Standards and Technology or Purdue 
University, nor does it indicate that the products are necessarily the best available for the purpose. 
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Figure 2-4. (a) Schematic view of LGCC experiment with adapted finite-difference nodes; (b) initial 
temperature of finite difference simulation, i.e., T (x,t = 0); and (c) temperature at the bottom of the LGCC 


experiment, i.e., T (x = 1,t), as a function of time. 
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Table 2-1. Thermal properties of the thermally conductive pad, foam and Pyroceram code 9606 


Material k (W/(m:K)) p (kg/m’) C (kI/(kg-K)) 
Thermal Pad 3.0 309 850 
Foam (Williams & 
Aldao, 1983) 0.03 20 1340 
Pyroceram Code 9606 
(Yaghoob Farnam, Bentz, á 2600 900 


Sakulich, et al., 2014) 


Note: *The thermal conductivity of Pyroceram code 9606 as a meter bar material was calculated as a function of 
temperature (Yaghoob Farnam, Bentz, Sakulich, et al., 2014). 


The first step in the numerical approach was to discretize the experimental setup by a finite 
difference method using an appropriate grid spacing size, Axof 1 (mm) and time step, At of 
0.05 (sec). The initial temperature of the entire experimental setup was set equal to the ambient 
temperature T (x, t=0) = 24 (°C), as displayed in Figure 2-4b. The temperature at the bottom of 
the LGCC experiment, T (x = / mm, t) (see Figure 2-4a), varied in the numerical simulation as a 
function of time according to the LGCC experimental protocol (Figure 2-4c). A heat convection 
coefficient Acom = 100 W/(m?:K) is employed to simulate the heat transfer between the air and the 
foam on the top (Incropera et al., 2007). Even though the insulating foam has quite a small thermal 
diffusivity parameter, significant heat energy 1s still transferred to the environment, resulting in a 
slight temperature differential between the top of the foam and the ambient environment. The 
relevant thermal properties of the thermally conductive pad, foam, and Pyroceram code 9606 used 
in the modeling are listed in Table 2-1. 

Figure 2-5 displays a flow chart of the one-dimensional explicit finite difference method 
adopted to simulate the saturated mortar specimen containing de-ionized water and NaCl solution 
with various concentrations. First, the thermal properties of components of the mortar specimen, 
temporal and spatial step sizes, and thermal initial and boundary conditions are determined. All of 
the discretized layers (1.e., x = 1 mm to x = 204 mm) are employed to simulate the heat transfer 
for the LGCC experiment; however, only the finite layers of the mortar specimen (1.e., x = 32 mm 
to x = 83 mm) are investigated in this figure. During phase transformation of the pore solution in 
the mortar specimen filled with de-ionized water, two approaches of consideration of either a 


discrete pore size distribution or a continuous pore size distribution were employed. The approach 
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of using a distribution of continuous pore sizes introduces a progressive ice formation/melting in 
the pores that are simulated using equations 2-3 and 2-4 and Figure 2-1; however, the volume 
fraction of pores with the size greater than 5nm (vr = 60%) is considered to simulate the 
instantaneous ice formation/melting occurring in the approach of using discrete pore sizes. 


Additionally, the progressive fraction of produced/melting ice is calculated using the lever rule. 
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Figure 2-5. Numerical algorithm of finite difference strategy using heat transfer equation to 
simulate the thermal behavior of a saturated mortar specimen. 
2.5 Results and discussion 


In this section, mortar specimens saturated with water at 75 %, 85 %, 95 %, and 100 % 


degrees of saturation (Ds) and the effect of pore size distribution are numerically investigated. 
2.5.1 Fully saturated mortar specimen 


Two numerical models, with either a continuous or discrete pore size distribution, are 


investigated in this section. Figure 2-6 shows the experimental and numerical results for the 


49 


thermal behavior of mortar specimens that were saturated (1.e., Ds = 100 %) with water. The heat 
loss coefficient, loss is assumed to be a constant value of 60 % in this figure. The model with a 


discrete pore size distribution only considers the instantaneous freezing of the pore solution that 


can freeze (vp =60%). In contrast, the model with the continuous pore size distribution also 


considers an additional process of gradual freezing of the pore solution (vp = 72%), as displayed 


in Figure 2-6(c). Figure 2-6(a) indicates the numerical and experimental temperature profile for the 
saturated mortar specimen at the bottom layer in the LGCC experiment setup (1.e., x = 32 mm). A 
nearly instantaneous temperature rise that occurred at the moment of freezing can be observed, 
demonstrating the instantaneous freezing due to the temperature of the first finite difference layer 
of the mortar specimen (1.e., x = 32 mm) reaching the freezing temperature of (undercooled) water 


in the mortar specimen (Tp =—6.1 °C). However, the model with a continuous pore size 


distribution considers that the remaining amount of unfrozen water in the pores transforms to ice 
gradually until the temperature reaches —35 °C. 
The model with a discrete pore size distribution considers the ice melting to occur by 


absorbing sufficient heat at the melting temperature of the large frozen pores observed in the 
experiment (7,, =0 °C). In the model with a continuous pore size distribution included, the ice 


formed in the pores is considered to melt gradually starting at —35°C (Z. Sun & Scherer, 2010b), 
based on the measured absorption isotherm for the mortar. 

The heat flow is obtained using the numerical simulation to evaluate the role of pore size 
distribution and compared to the experimental data shown in Figure 2-6(b). The formation of ice 
in the pore solution results in an exothermic peak, which is representative of the latent heat released 
during a freezing cycle. In the model with a discrete pore size distribution, the exothermic peak is 
considered to occur at Ts. Subsequently, it ceases when the entire amount of latent heat has been 
emitted to the surroundings. The endothermic peak begins as a gradual process at 0 °C until all of 
the previously formed ice melts inside the frozen pores. 

Conversely, the exothermic peak is extended to the end of the freezing cycle (—35 °C) due 
to gradual ice nucleation inside the smaller pores in the model with a continuous pore size 
distribution. For the case of melting, the endothermic peak is considered to occur gradually as a 


function of temperature and the pore size. 
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Figure 2-6. Experimental and numerical results for mortar specimen fully saturated with water (a) 
temperature profiles at x = 32 mm (see Figure 2-1 for the definition of x); (b) predicted heat flow; 
(c) volume of pore solution undergoes a phase transformation at the bottom surface of the mortar 


specimen. 


Therefore, the melting curve extends progressively to 0°C, owing to the broad range of 
pore sizes in the model with a continuous pore size distribution. It is concluded that the 
consideration of pore size distribution can reasonably be neglected during the freezing process due 
to undercooling. In contrast, the melting of formed ice indicated a gradual process as the 


temperature increases in both the experimental data and the model with a continuous pore size 


distribution. 
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2.5.1.1 Partially saturated mortar specimen 


The amount of heat released during freezing (AH A | was obtained using the numerical 


simulation for mortar specimens saturated at different degrees of saturation with the discrete pore 
size distribution model and is compared with experimental results in Figure 2-7. For partially 
saturated specimens, the value of vr was calculated using vr=Ds-vw, where vw1s the volume fraction 
of pores with a size of less than 5 nm (vw = 40 %). The results for these partially saturated mortar 
specimens are the experimental investigation of this work, whereas the experimental data of the 
fully saturated mortar specimen was already published (Yaghoob Farnam, Bentz, Sakulich, et al., 
2014). Two heat dissipation coefficients of 40 % and 60 % are considered. The coefficient of 
variation for the experimental results is obtained at 8.6 %. The numerical simulation predicts 


greater heat release than that obtained in the experiment due to further experimental imperfections. 
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Figure 2-7. The amount of heat released during freezing (AH;,, ) for the mortar specimen 
saturated at different degrees of saturation (Ds). 


The amount of heat released during freezing (AH,,,) was obtained using the numerical 
simulation for mortar specimens saturated at different degrees of saturation with the discrete pore 
size distribution model and is compared with experimental results in Figure 2-7. For partially 
saturated specimens, the value of vr was calculated using vr=Ds-v,, where vy is the volume fraction 


of pores with a size of less than 5 nm (vw = 40 %). The results for these partially saturated mortar 
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specimens are the experimental investigation of this work, whereas the experimental data of the 
fully saturated mortar specimen was already published (Yaghoob Farnam, Bentz, Sakulich, et al., 
2014). Two heat dissipation coefficients of 40 % and 60 % are considered. The coefficient of 
variation for the experimental results is obtained at 8.6 %. The numerical simulation predicts 


greater heat release than that obtained in the experiment due to further experimental imperfections. 
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Figure 2-8 The amount of heat released during freezing (AH -») for the mortar specimen saturated 
at different degrees of saturation (Ds). 


2.6 Summary and conclusion 


In this paper, a one-dimensional finite difference numerical model was used to predict the 
macroscopic freeze-thaw behavior of air-entrained mortar specimens. The effective thermal 
properties of the composite mortar were estimated using homogenization techniques. The role of 
curvature, owing to a broad range of pore sizes, was considered in calculating the volume fraction 
of freezable pore solution exposed to freezing/thawing cycles using measured absorption- 
desorption isotherms. It was concluded that the role of pore size (or curvature) on the macroscopic 
behavior of the air-entrained mortar specimen was negligible during freezing due to the quantity 
of larger pore sizes in realistic mixtures and undercooling. In contrast, the role of curvature had a 
considerable impact on the macroscopic behavior of the frozen mortar specimen during melting. 


The lever rule approach, derived from a phase diagram of the NaCl-water solution, and 
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undercooling were adopted in the numerical model. It was concluded that this model can simulate 
the freezing and thawing process of mortar specimens saturated with water or various NaCl 
solutions to predict the thermal behavior of mortar specimens at various degrees of saturation or 
saturated with various concentrations of NaCl solutions. 

The computational results were compared to the experimental ones obtained for mortar 
specimens saturated with NaCl solution using a low-temperature longitudinal guarded comparative 
calorimeter (LGCC). A lower amount of heat release (or freezable fraction of pore solutions) was 
observed in the experiment than the theoretical value predicted based on the measured desorption 
isotherm. The difference may be mainly due to experimental conditions allowing significant heat 
dissipation within the LGCC experiment. To justify the experimental under-estimation of heat 
release, two heat loss coefficients of 40 % and 60 % were evaluated to validate the numerical 
results. Accordingly, a better agreement was exhibited between the numerical results and the 


experimental data. 
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3. ENHANCING FREEZE-THAW PERFORMANCE OF 
CEMENTITIOUS COMPOSITES THROUGH INCORPORATION OF 
PHASE CHANGE MATERIAL (PCM) 


This chapter contains work that was originally published in Elsevier as “Hadi S. Esmaeeli, 
Yaghoob Farnam, John E. Haddock, Pablo D. Zavattier1, W. Jason Weiss. Numerical analysis of 
the freeze-thaw performance of cementitious composites that contain phase change material 
(PCM). In Journal of Materials and Design. Volume 145., pp. 74-87. Elsevier, 2018.” The original 
article has been used with permission as stated below. 

Reprinted by permission from Elsevier, In: Journal of Materials and Design. “Numerical 
analysis of the freeze-thaw performance of cementitious composites that contain phase change 
material (PCM).” Hadi S. Esmaeeli, Yaghoob Farnam, John E. Haddock, Pablo D. Zavattieri, W. 
Jason Weiss, COPYRIGHT (2018). 


3.1 Introduction 


Two main factors affect the freezing and melting behavior of PCM impregnated within LWA: 
(1) the rate of heat release or absorption of bulk PCM, and (2) pore size distribution inside LWA. 
Low-temperature differential scanning calorimeter (LT-DSC) is a standard method to analyze the 
thermal energy storage capacity of PCM and to determine the thermal properties of PCM, 
including the transition temperature range and 4H;-temperature relationship (Castellon et al., 2008). 
According to the Gibbs-Thomson equation, the pore size causes a capillary pressure in the solution 
inside the pore, which leads to depressing the freezing temperature of PCM, depending on the pore 
size (Z. Sun & Scherer, 2010b). Therefore, the size of pores in the LWA can alter the freezing 
temperature of PCM and the corresponding 4H;-temperature relationship. In this work, we employ 
paraffin oil as the PCM. However, it is worth mentioning that undercooling does not occur to a 
significant extent for this particular PCM, as paraffin is chemically stable with a low vapor pressure 
(C. Liu et al., 2010). 

Improving the freeze-thaw performance of concrete pavements would have significant 
economic impacts. The American Society of Civil Engineers issued a grade of “D” to America’s 
road infrastructure, as some $67 billion was spent annually on the repair of damaged roads in 2013 
(Herrmann, 2013). It 1s believed that PCM incorporation can be used as an effective method to 
increase the thermal energy storage capacity by reducing the depth and time of freezing that a 
concrete pavement can experience, thus minimizing concrete pavement damage. A finite- 


difference approach is implemented to simulate the concrete pavement considering realistic 
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thermal (environmental) boundary conditions. PCMs with three transition temperatures, which are 
all designed above the freezing point of pore solution, are considered. Year-long simulations 
throughout the United States are performed to demonstrate the impacts of transition temperatures 
and AH;-temperature relationship of PCMs on reducing the time and depth of freezing within the 
pavements. The outcomes provide new insights that are needed to ascertain the optimal PCM 


characteristics that mitigate the risk of freeze-thaw cycling in concrete pavements. 


3.2 Numerical modeling 


This study aims to determine the thermal behavior of cementitious composites such as mortar 
and concrete specimens in which the incorporated PCM undergoes phase transition during cooling 
and heating cycles. The heat release/absorption is determined as well as the corresponding fraction 
of PCM that can freeze/melt within the LWA pores. Consequently, the temperature profile and 
associated heat flow within the specimen can be calculated by solving the heat transfer equation 


shown in equation 3-1 (Hadi S. Esmaeeli et al., 2017). 


€°T (x, t) 
ax? 


OT (x,t) 


. ; ; 3-1 
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Ke (Tr). 
where T (x,t) is the temperature of the specimen at the location of x (mm) and time of t (sec) (°C), 


k, is the thermal conductivity of specimen at temperature T[W/m.K)], pc is the density of the 


specimen at temperature T (kg/m’), C 4 is the heat capacity of the specimen at temperature 


T [JXke.K)], EE is the rate of heat release or absorption associated with freezing or melting 


of PCM, as described in equation 3-2a [JAm?.sec)], a taren 1o, CODSIders the rate of heat dissipation 
during phase transition of PCM (to the environment), as described in equation 3-2b [JAm?.sec)], 
and Foon, is the rate of heat dissipation throughout the experiment due to the temperature gradient 


between the specimen and surroundings (or vice versa) through convection as described in 


equation 3-2c [JAm’.sec)]. The heat sink term representing the heat dissipation during phase 


transition of PCM is a fraction of heat release/absorption which is 


7 atent-loss — flaten- gen’ hiren toss 
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determined using the dissipation coefficient of /jatent-loss (<1) obtained from experiments (Hadi S. 
Esmaeeli et al., 2017; Yaghoob Farnam et al., 2017; Yaghoob Farnam, Krafcik, et al., 2015). 
Equation 3-2c describes a convection heat term that represents the heat dissipation to the 
surrounding using the convection coefficient, Acony-diss which is also obtained empirically (Yaghoob 


Farnam et al., 2017; Yaghoob Farnam, Krafcik, et al., 2015). 











| Ovr (T) 
d latent-gen E AH + ' P Pcm ` VYPCM `’ za 3-24 
| Ovr (T) 
d latent-loss ~ AH Po ree ee ' P Pcm ` VPCM ` 7 A 3-2b 
r -h éT (x,t) 

conv conv —diss ox 3-2c 


where ppcm 1s the density of phase change materials (PCM or pore solution) within the specimen 
in this study (PCM or pore solution) (kg/m), vecuis the volume fraction of phase change materials 
(PCM or pore solution), vr is the volume fraction of phase change materials (PCM or pore solution) 
that can freeze or melt at temperature (T). The variable vris calculated by normalizing the amount 
of heat release/absorption at the temperature T respect to the total amount of heat that PCM can 
release or absorb. 

In this study, two phase change materials, (1) water solution absorbed within the pore 
structure of mortar and concrete specimens, and (2) PCM absorbed within the pore structure of 
LWAs, each with distinct thermal properties, are studied in two small-scale longitudinal guarded 
comparative calorimeter (LGCC) and large-scale (slab) tests. PCM is designed to experience phase 
transition at a temperature above the freezing point of pore solution, where the freezing of pore 
solution depends on the thermal (environmental) and saturation conditions of specimens. An 
LGCC test was conducted on a mortar specimen, which was prepared using LWAs fully saturated 
with PCM and then maintained in a chamber with a controlled relative humidity of 50 + 1 % 
(Yaghoob Farnam, Krafcik, et al., 2015). A slab test was performed on a concrete specimen, which 
was prepared using LWAs fully saturated with PCM, then moist cured for 14 days and air-cured 
until the testing time (Yaghoob Farnam et al., 2017). Under these conditions, the internal relative 


humidity of the slab was reduced to 81 % (Jun Zhang et al., 2012). Esmaeeli et al. (Hadi S. 
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Esmaeeli et al., 2017) predicted that the impact of the solution inside the pores smaller than 5 nm 
(corresponding to a relative humidity less than 80 %) on freezing was relatively insignificant. It 
can be concluded that the pore solution inside the mortar and concrete is less likely to experience 
phase transition, thus considering the PCM as the only phase change material in this study. This 
prediction is consistent with the experimental results of the LGCC and slab tests where no phase 
transition events associated with freezing of pore solution were observed (Y. Farnam et al., 2015; 


Yaghoob Farnam et al., 2017). 


3.3 Numerical methodology and applied boundary conditions 


An in-house one-dimensional explicit finite difference code has been developed using a C 
programming language to predict the thermal response of cementitious composites with and 
without freeze-thaw cycling applied to the PCM. Depending on the test (e.g., LGCC or slab test), 
the experimental setup is discretized spatially and temporally using a convergent grid spacing size, 


Ax, and time step, At. The cementitious composites had effective thermal properties k-(T), p-(T), 


Ce r) and experienced a local rate of heat release/absorption per unit volume within each finite- 


difference grid . In the LGCC test, the bottom surface is exposed to heat conduction with 


d latent—gen 
a heat sink term. In the slab test, the top surface of the pavement was exposed to convection with 
the ambient temperature with convective coefficient Aconv-aiss. At time t = 0 sec, the temperature in 
the cementitious composites is considered to be uniform and equal to ambient temperature, Tamb, 


as described in equation 3-3. 
T (x,t =0)=Toyn, 3-3 


Further details on the simulation methodology can be found in (Hadi S. Esmaeeli et al., 2017). 


3.4 Freezable fraction of PCM, vr (T) 


Paraffin oil (Ci4-Cie) is used as a commercially available PCM in this study, which is a 
petroleum-based material with a specific density of 0.77 and a vapor pressure less than 0.01 mmHg. 
Figure 3-1 shows the thermal energy storage analysis to characterize the thermal properties of bulk 
PCM, including transition temperature range and the associated AH. A TA instrument Q2000 LT- 


DSC instrument with temperature range varied between -40 °C to +24 °C with a cooling rate of 
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+2 (°C/h) and heating rate of +4 (°C/h) (identical to the LGCC experiment protocol) is used. The 
paraffin oil specimens with a mass of 10 + 2 mg are weighed into a T-zero stainless-steel high- 
volume pan (with 100 uL volume capacity) and sealed hermetically. Dry nitrogen is used as the 
specimen purge gas and cooling gas with the flow rate regulated at 50 (mL.min’). Figure 3-1(a) 
shows the DSC thermogram of PCM, where the heat flow of PCM is plotted as a function of 
temperature. In the DSC thermogram, the main exothermic peak labeled as “A” represents solid- 
liquid phase transition associated with freezing of PCM at the temperature of +4 °C. In comparison, 
minor peak labeled as “B” corresponds to the solid-solid phase transition of PCM at a temperature 
of -27.5 °C. 

On the other hand, minor endothermic peak labeled as “C” represents the phase transition 
associated with solid-solid of PCM at a temperature of -28.4 °C. In contrast, major peak labeled as 
“D” corresponds to solid-liquid phase transition associated with the melting of PCM at 
temperatures of +5.1 °C. The heat released during the cooling cycle and heat absorption during the 
heating cycle are calculated and shown in Figure 3-1(b). For this study, only the phase transition 
of PCM at the “A” and “D” peaks are considered since (1) a PCM with an initial freezing 
temperature somewhat above the initial transition temperature of pore solution (assumed O °C in 
this study) has good performance in ice and snow melting/removal applications (Yaghoob Farnam 
et al., 2017; Yaghoob Farnam, Krafcik, et al., 2015) and also leads to improvement in freeze-thaw 
performance of cementitious materials (Dale P. Bentz & Turpin, 2007; T.-C. Ling & Poon, 2013; 
A.R. Sakulich & Bentz, 2012); and (2) AHş measured at these two peaks are significantly higher 


comparing to 4H; measured at the “B” and “C” peaks. Hence, a modified version of fusion of 


enthalpy “ AH A ” is defined as a function of temperature for freezing and melting processes, as 
shown in Figure 3-1(b). For the freezing process, it is observed that AH F begins to release once 


temperature reduces to +4°C (1.e., the initial freezing temperature of PCM, i) and increases 


gradually until temperature reduces, to approximately -10 °C (.e., the final freezing temperature 


of PCM, T! ). As the temperature drops below -10 °C, no considerable change in AH F is observed. 


For the melting process, no change in AH F is detected until temperature increases to -8 °C (1.e., 


the initial melting temperature of PCM, T ). At this temperature, AH F begins to retrieve the 
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released heat gradually until the total heat is absorbed at a temperature of +5.1 °C (.e., the final 


melting temperature of PCM, Ti ). On the other hand, the manufacturer of this particular PCM has 


reported a freezing temperature of 1.7 °C and a melting temperature of 4.5 °C. It is necessary to 
note that the difference between the experimental and reported transition temperatures may be 
caused as the DSC results highly depend on the cooling and heating rates (Watson et al., 1988). 


Although the slower cooling and heating rates more evenly affect the specimen which can lead to 


reduce the calculated i and Ti to the reported freezing and melting temperatures. The cooling 


and heating rates employed in the LT-DSC test is consistent with the LGCC experiment protocol 
which enables us to validate the numerical model by comparing the numerical results with 


experimental data obtained for mortar specimen containing PCM subjected to freeze-thaw cycles. 
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Figure 3-1. (a) Thermal storage analysis of PCM using an LT-DSC device, (b) the amount of 
released and absorbed heat, 4H;, calculated at exothermic peaks of “A-B” and endothermic 
peaks of “C-D” using LT-DSC result, and the amount of released and absorbed modified heat, 


AH f , used in the numerical analysis. 


The pore size distribution of LWA is also a primary factor that can impact the freezing and 
melting processes of PCM. The LWAs contain a pore structure with a broad range of sizes, which 


can alter the freezing temperature of liquid PCM (Yaghoob Farnam, Krafcik, et al., 2015). The 
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Gibbs-Thomson equation can be used to correlate the freezing temperatures of a liquid PCM inside 


an LWA pore to the pore size, as shown in equation 3-4. 


YL/s * UL 


Tr = IM — 21M Nor 3-4 


where y;,/s1s the energy of the liquid/solid PCM interface (~ 0.0488 (J/m’)) (Davies, 2012), r is 


the pore radius (size), vz is the molar volume of the liquid PCM (= 0.0002761 (m7/mol)) (Huggins, 
1941), Ahysis the molar enthalpy of melting of PCM (~ 3850 (J/mol) (Jackson & McKenna, 1990), 


Tr is the freezing temperature of PCM inside pore size r, and Ti is the final melting temperature 


of an infinitely large solid PCM (Ti is obtained from LT-DSC analysis) (Parks & Todd, 1929). 


Figure 3-2(a) displays the correlation between the pore size and the temperature at which PCM 
inside the LWA pore freezes, 7; (r). Once the temperature reaches the associated freezing 
temperature, solid PCM begins to form inside the pore with radii of r. This process is reversed 
during melting. 

Figure 3-2(b) illustrates the relationship between vr and r obtained from the empirical 
measurement of the size distribution of LWA pores (Castro, Keiser, et al., 2011). It is observed 


that 84 % of the total PCM by volume absorbed in the LWA pores with sizes larger than 29 nm 


can freeze once the temperature drops to +4 °C (Ti). In comparison, 99 % of liquid PCM by 


volume absorbed in the LWA pores with sizes larger than 2.13 nm can freeze as the temperature 


drops further to -10 °C (TZ ). Consequently, the volume of freezable PCM determined without 


considering the effect of the size distribution of LWA pores is underestimated by about 15 %, 
compared to the approach where the size distribution effect of LWA pores is considered. Esmaeeli 
et al. observed that the approximate 16 % underestimation of the volume of freezable pore solution 
absorbed inside the mortar has a negligible impact on the macroscopic freezing response of 
specimen (Hadi S. Esmaeeli et al., 2017). Consequently, the work reported in this paper considers 
the pore size of LWA to have a negligible contribution to the macroscopic thermal behavior of 
cementitious materials containing PCM. 

To calculate vr, two approaches are evaluated in this study: (1) a model considering the 


instantaneous phase transition of PCM and (2) a model considering a gradual phase transition of 


6l 


PCM. The first approach only considers the instantaneous freezing or melting of PCM at 


temperatures of I or Ty, , respectively, whereas the second approach considers a gradual process 
of phase transition of PCM. In a second approach, the freezing process is considered to initiate L 
and gradually continues until the total heat is released Ti . In contrast, the melting process is 


considered to initiate l; and continues gradually until the total released heat is absorbed at Ti ; 


as shown in Figure 3-1(b). 
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Figure 3-2. (a) The effect of pore size (r) on the freezing temperature of liquid PCM using the 
Gibbs-Thomson equation, (b) volume fraction of liquid PCM that can freeze (vr) as a function 
of pore size (7). 


3.5 Homogenization of the thermal properties of cementitious composites 


This section discusses how the thermal properties of a mortar or concrete specimen 
composed of three constituents with distinct thermal properties (1.e., dry mortar or concrete, water 
solution absorbed within the pores of the specimen, and PCM incorporated within the LWA pores) 
can be determined using homogenization techniques such that they can be used in the finite 
difference approach. 

Depending on the saturation state and temperature of the specimen, pores may contain 


various constituents, including air, pore solution, liquid PCM, or solid PCM. The amount of 
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solution and air in the specimen depends on its degree of saturation (Ds), where the amount of 
liquid and solid PCM depends on the temperature of the specimen. Two PCM dosages of 
145 (kg/m’) and 170 (kg/m”) are used to make mortar and concrete specimens tested in the LGCC 
and slab tests, respectively. Subsequently, the change in the amount of air, pore solution, liquid 
PCM, or solid PCM can substantially alter the effective thermal properties of the specimen due to 


considerable differences between the thermal properties of constituents, as outlined in Table 3-1. 


3.5.1 Thermal conductivity, ke 


The rule of the mixture is used to determine the effective thermal properties of pores, kp, 
as a function of the volume fractions of each constituent (1.e., Kair, kw, kip, and ksp), which is 
outlined in equation 3-6. Afterward, the effective medium theory (EMT) can be used to determine 
the thermal conductivity of a cementitious composite, ke, composed of the dry cementitious 
composite as matrix and pores as inclusions of composite and using equation 3-5 (Levy & Stroud, 


1997). 


air* air 


P P P 
AR E a a + kiplveyg—Ve + K pV rq 3.5 


21 —v? — veya) Kary + (L+ 2v? + 2vP 4) kp 3-6 


k.=k 
2+v? +Viy,}k 


c T “dry” P P 
dT l—vo —Viwa) Kp 


where kary is the thermal conductivity of dry cementitious material [W/Am.K)], ye is the volume 
fraction of pores in the cementitious material respect to the total volume of pores (vi pg am ), 
aon is the volume fraction of pores in the LWAs respect to the total volume of pores 


(vi + a ), kair is the thermal conductivity of air [WAm.K)], Vair=1-vw is the volume fraction of 


pores of cementitious material filled with air respect to the total volume of pores in the 


cementitious material (v), kw 1s the thermal conductivity of water solution [WXm.K)], vw is the 


volume fraction of pores of cementitious material filled with water solution concerning the total 
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volume of pores in the cementitious material (vi ), kıp is the thermal conductivity of liquid 


paraffin [W/m.K)], and ksr is the thermal conductivity of solid paraffin [W/m.K)]. 


Table 3-1. Thermal properties of constituents of mortar and concrete composites containing 
water solution and PCM 


Material k (W/m.K) p (kg/m’) œ [KJ/kg.K] 
Air 0.023 1.35 1.005 
Dry mortar 1.78 + 0.3 2162 + 150 950 + 100 
Dry concrete 1.70.1 2200 100 880 + 70 
Water solution 0.52 + 0.018 998 +5 4183 +5 
Liquid paraffin 3.1 +0.7 766 + 340 2981 + 190 
Solid paraffin 32 0.7 865 + 148 2604 + 248 


3.5.2 Density, pc, and specific heat capacity, C A 


The law of mixtures is employed to predict the effective density and heat capacity of the 


cementitious composite, pc and C A . The effective density is estimated using equation 3-7. 


P, = Pay (l-Vp)+ Pp “Vp 3-7 


where paris the density of dry cementitious composite (kg/m”’), pp is the density of the materials 


inside the pores (i.e., pore solution and PCM) (kg/m”), which is calculated using equation 3-8. 


= P P P 
Pa~ Pir Vae Ve TPN Ve Poy, — Vp )+ Psp YF 3-8 


where pairis the effective density of air (kg/m’), pw is the density of water solution (kg/m”°), pip is 
the density of liquid paraffin (kg/m”), and psp is the density of solid paraffin (kg/m”). 

The process of determining an effective heat capacity is conceptually similar to 
determining an effective density (Hadi S. Esmaeeli et al., 2017; L.-P. Zhou et al., 2010). The 


effective heat capacity for a cementitious composite is calculated using equation 3-9. 
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C? - p, =Chy* Pay (vp }+ CB: Pp vp ae 


where C a is the heat capacity of the dry composite [Jkg.K)] and C P is the effective heat 


capacity of the materials inside the pores, where C A is calculated using equation 3-10. 
CP pp =C pu, Van Vi +C. p, v evi +CP b? - +c? 
P ` Pp = gir-Pair ` Vair “Vo T “w-Pw' Vw" Ve + ULP-PLp-Wiwa ~ VF }* Usp-Psp-YF 3-10 


where C is the heat capacity of air [JXkg.K)], CE is the heat capacity of water 


solution [JXkg.K)], Cj» is the heat capacity of liquid paraffin [JXkg.K)], and C¿p is the heat 
capacity of solid paraffin [Jkg.K)]. 


3.6 Simulation results 
3.6.1 Sensitivity analysis 


Table 3-1 outlined the input variables representing the thermal properties of individual 
constituents of cementitious composites that contain PCM. The range of input variables indicates 
an uncertainty that may influence the thermal behavior of the specimen as the output variable. 
Hence, a sensitivity analysis 1s required to investigate how the uncertainty in the output variable 
can be apportioned to different sources of uncertainty in the input variables. In this analysis, the 
numerical analysis of the LGCC test is used to study the sensitivity of the input variables by 
employing the minimum and maximum values for input variables to quantify the heat flow 
response of a mortar (1.e., output variable). Following the experimental work by Farnam et al. (Y. 
Farnam et al., 2015), the LGCC test was designed to quantify the temperature profile and 
corresponding heat flow within a mortar specimen. More information on the mixture proportioning, 
specimen preparation, and conditioning of the mortar can be found in (Y. Farnam et al., 2015). 
Figure 3-3(a) displays the LGCC test setup, where a one-dimensional heat transfer is conducted 
by providing a heat sink at the bottom, longitudinal thermal insulation on the sides and the top, 


and foam insulation around the specimen to minimize the specimen/environment heat transfer. 
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The first step in the numerical approach is to employ the finite difference method using an 
appropriate grid spacing size, Ax of 1 mm, and time step, At of 0.05 sec. Figure 3-3(b) shows the 
temperature at the bottom of the LGCC test, T (x = 1 mm, t) that imposes one thermal cycle by 
varying the temperature as a function of time within the range of +24 °C to -35 °C, according to 
the LGCC test protocol. At x = 2 mm, the temperature in the first finite-difference grid of pad is 
determined by imposing a heat flux boundary condition to include the effect of conductive heat 
transfer, as described in equation 3-11. 


a°T(x,t) B T 


Kpad > 2 = Ppad C paa T 


The heat convection coefficient at the top surface (x = 204 mm) is considered to be 
approximately Aconv-aiss % 100 [WAm7.K)] to simulate the heat transfer between the air and the foam 
on the top (Incropera et al., 2007). At x = 203 mm, the temperature in the top finite-difference grid 
of foam is determined by imposing a heat flux boundary condition to include the effect of 


convective heat transfer occurring between foam and environment, as described in equation 3-12. 


OT (x, t) 7 


E Kain ET B Me oes l Tomt ()- T(x, r )| 3-12 


A considerable temperature difference between the specimen and surroundings throughout 
the experiment, where a heat flux boundary condition can be taken into account through the lateral 


side of the experimental setup, as described in equation 3-13. 


e) hy aos Fole) - Te) 3-13 


Additionally, a difference between the measured heat released in the experiment and the 
associated AH; of incorporated PCMs was observed (Hadi S. Esmaeeli et al., 2017; Yaghoob 
Farnam, Bentz, Sakulich, et al., 2014; Yaghoob Farnam, Krafcik, et al., 2015). Volume fractions 


of PCM (vim), pore solution (vw), and air (Vair) concerning the total volume of the pores 


(V F vy ) are calculated as 42 %, 23 %, and 35 % in the LGCC test, respectively (Yaghoob 


Farnam, Krafcik, et al., 2015). 
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Figure 3-3. (a) Schematic presentation of the LGCC experiment with imposed finite-difference 
nodes at certain levels, (b) temperature profile at the bottom of LGCC test, 1.e., T (x = 1 mm, t), 
adapted from (Hadi S. Esmaeeli et al., 2017). 


Figure 3-4 shows the heat flow response of mortar in the LGCC test, where eight 
parameters Sc7, Sc2, Sai, SH2, Leri, Ler, Heri, and Hpr2 are introduced to characterize heat flow. 
For the cooling cycle, two parameters Sc; and Sc2 control the slope of the curve before initiation 
and after completion of the freezing process of PCM. In comparison, two parameters Hpr; and Lpr; 
control the maximum heat flow and exothermic peak time during the freezing process of PCM. 
For the heating cycle, the two parameters Sq; and Sp2 control the slope of curve initiation and after 
completion of the melting process of PCM. In contrast, the parameters Hpr2 and Lpr2 control the 


maximum heat flow and exothermic peak time during the melting process of PCM. 
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Figure 3-4. The categorization of heat flow response of mortar into six parts of (C/) cooling 
before phase transition of PCM, (PT1) cooling during phase transition of PCM, (C2) cooling 
after a phase transition of PCM, (H1) heating before phase transition of PCM, (P72) heating 
during phase transition of PCM, and (H2) heating after a phase transition of PCM. 


Figure 3-5 shows the results for sensitivity analysis of seven input variables, kary, Pary, C A : 


kip, ksp, Miatent-los, and Acony-diss. For simplicity, this analysis considers the approach of the 
instantaneous phase transition of PCM to simulate the freezing and melting processes. It is also 
necessary to mention that the water transport occurring during the freeze-thaw cycling is neglected 
as no phase transition associated with ice formation inside the pores is expected (G. W. Scherer & 
Valenza, 2005). 

The coefficient Aconv-aiss can describe the heat energy transferred between specimens and 
the environment through convection (Incropera et al., 2007). The typical convection coefficient 
between air and objects is considered in the range of 1-25 [W/(m7.K)] [37, 38]. As two layers of 
foam and longitudinal insulators are placed on the sides of the LGCC apparatus to minimize the 
heat transfer to the environment, hconv-aiss can be considered much less than 1 [W/(m7’.K)]. The 
increase Of Acony-diss WOuld increase the temperature gradient within the specimen, thereby 
increasing the amount of heat flow. Figure 3-5(a) also shows that a slight change in Aconv-adiss from 


0 to 1 [W/(m7.K)] significantly increases the four slope parameters Sc7, Sco, Sui, and Sp2. 
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The variation and the range of Kary pary and C A y are respectively shown in Figure 3-5(b), 


Figure 3-5(c), and Figure 3-5(d). The increase of kary would either (1) decrease the temperature 
gradient within the specimen, which may decrease the heat flow, or (2) increase the effective 
thermal conductivity of specimen, which may increase the heat flow. The sensitivity analysis result 
shows that the increase of kary from 1.5 to 2.1 [W/(m.K)] results in an increase of heat flow, and 
therefore the slope parameters Sc7, Sco, Sai, and Sx2 are decreased; this change has a relatively 
smaller impact than does the variation of /Acony-aiss. Furthermore, parameters Hpr; and Hp72 decrease, 


whereas parameters Lpr; and Lp72 increase. Conversely, the increase of pary from 2012 to 


2312 (kg/m?) and C? from 850 to 1050 [JAkg.K)] would result in a decrease of a temperature 
dry 


gradient within the specimen and therefore decrease the heat flow throughout the experiment 


except during the phase transition of PCM. As such, the increase of pary and C, 5 tend to decrease 


the slope parameters Sc7, Sc2, SH1, and Sp2. 
Figure 3-5(e) and Figure 3-5(f) show how the variation of thermal conductivities of liquid 


and solid PCM (krr and ksp) can influence heat flow. The density and heat capacity of liquid and 


solid PCM (prp,C;p psr Cip) have a negligible influence on heat flow due to a relatively low 


volume of cementitious composites occupied by PCM. The sensitivity analysis results show that 
an increase of kzp from 0.21 to 0.30 [W/(m.K)] increases the heat flow before initiation of freezing 
or after completion of the melting of PCM. On the other hand, an increase of ksp from 0.21 to 
0.30 [W/(m.K)] increases the heat flow after completion of freezing or before initiation of the 
melting of PCM. Therefore, the increase of kzp only decreases two slope parameters Sc; and Sz, 
while the increase of ksp only decreases two slope parameters Sc2 and Sz7. 

Figure 3-5(g) shows how the input variable /jatent-ioss associated with the heat energy 
dissipation during the phase transition of PCM significantly influences the heat flow. The 
sensitivity analysis result indicates that the variation of /jatent-ioss from 40 % to 60 % would result 
in a significant decrease of the area under the exothermic and endothermic curves, and therefore 
decrease the slope parameters Hp7;, Lepr, Hp72, and Lp72. Table 3-2 summarizes the list of effective 


input variables indicated by sensitivity analyses on the corresponding heat flow parameters. 
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Figure 3-5. Assessment of the uncertainty of the input variables (a) convection dissipation 
coefficient, (b) thermal conductivity of dry mortar, (c) density of dry mortar, (d) heat capacity of 
dry mortar, (e) thermal conductivity of liquid PCM, (f) thermal conductivity of solid PCM, and (g) 
dissipation coefficient during phase transition of PCM, on the output variable (e.g., heat flow) of 
mortar specimen in the LGCC test. 
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Figure 3-5 continued 
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Table 3-2. The influencing input variables on the heat flow parameters of mortar composite 
containing PCM in the LGCC test. 


Parameters Sci Sc2 Apri Lpri SHI SH2 Hpr2 Lpr2 
Input hconv-diss  Nconv-diss kary kary hconv-diss — Nconv-diss kary kary 
variables kary kary Niatent-loss  Miatent-loss kary kary Niatent-loss  Miatent-loss 
Pdry Pdry Padry Padry 
ch œ, ch œ 
Kip ksp ksp Kip 


3.6.2 Small-scale LGCC test 


The numerical results for heat flow response of a mortar specimen containing PCM 
(incorporated within LWA) with dimensions of 25.4 mm x 25.4 mm x 50.8 mm tested in the LGCC 
test are compared to experimental data, as shown Figure 3-6(a). Two distinct numerical approaches, 
a model considering of instantaneous PCM phase transition and a model considering gradual PCM 
phase transition, are employed to investigate the impact of the transition temperature range of PCM 
on heat flow response. Table 3-3 shows the values for the thermal properties of mortar specimen 
constituents obtained from sensitivity analyses. The two dissipation coefficients Aconv-diss and Niatent- 
loss are also adjusted to be 0.7 [WAm’.K)], and 50 % for the LGCC test, respectively. 

Figure 3-6 (b) shows the freezing process of PCM results in an exothermic peak due to 
significant heat release during a thermal cycle. In the approach with instantaneous PCM freezing, 


the exothermic peak is considered to occur at T; and subsequently ceases when the total amount 


of heat has been released to the surroundings. Conversely, the approach with gradual PCM freezing 


considers the exothermic peak to be extended TZ due to a gradual process of solid PCM 


formation, as shown in Figure 3-6(b). For the case of instantaneous PCM melting, the endothermic 


peak is considered to occur at T, and subsequently ceases when the total amount of heat has 


been absorbed from the surroundings. In contrast, the endothermic peak is considered to be 


extended to T, as the previously formed solid PCM melts, as shown in Figure 3-6(c). The 


promising agreement between the experimental data and numerical results illustrates that 
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consideration of gradual freezing or melting of PCM over a narrow temperature range can 
reasonably simulate the thermal response of mortar specimens containing PCM exposed to one 


thermal cycle. 


3.6.3 Large-scale slab experiment 


A large-scale test with dimensions of 350 mm x 350 mm x 215 mm is simulated to assess 
the freeze-thaw performance of concrete slabs with and without PCM. A gradual phase transition 
of the PCM over a narrow temperature range is considered (see Figure 3-1(b)). The materials and 
mixture proportioning for making the concrete slabs are explained in Appendix 1. More 
information on the mixture proportioning, specimen preparation, and conditioning of the concrete 
slabs can be found in (Yaghoob Farnam et al., 2017). Figure 3-7(a) displays the slab test setup, 
where a one-dimensional heat transfer is conducted by providing a heat sink at the top, thermal 
insulation foam and wood frame with relatively low thermal conductivities of 0.03 [WAm.K)] and 
0.12 [WAm.K)], respectively, placed on the sides and bottom of the specimen (Sears et al., 1982; 
Williams & Aldao, 1983). The specimen domain is discretized spatially and temporally using layer 
size and a time step of | mm and 0.1 sec to obtain convergence in results, respectively. Therefore, 
this study considers 215 and 25 layers for the concrete specimen and base thermal insulator, 
respectively. Figure 3-7(b) shows the temperature at the top of the slab test, T (x = / mm, t) that 


imposes one thermal cycle by varying the temperature as a function of time within the range of 


+10 °C to -10 °C, according to the slab test protocol. Volume fractions of PCM ca ), water 


solution (vw), and air (Vair) concerning the total volume of the pores (ve T Ve ) are estimated to 


be 46 %, 45 %, and 9 % in the slab test, respectively (Yaghoob Farnam et al., 2017). 


Table 3-3. Thermal properties of constituents of mortar specimen in the LGCC test 


Material k (W/m-K) p (kg/m’) C” [kI/kg-K] 
Dry mortar 1.92 2080 870 

Liquid paraffin 0.21 670 3071 

Solid paraffin 0.21 800 2461 


73 


Heat Flow (W) 


0.01 





-20 -10 0 10 
T (x=32 mm, t} (°C) 


Heat Flow (W) 


Exp 
0.5 Num w/ instantaneous transition 
(W) —-—-—-—: Num w/ gradual transition 





-40 -30 -20 -10 0 10 20 30 
T (x=32 mm, t) (°C) 


(a) 


Heat Flow (W) 





0.01 


0 10 
T (x=32 mm, t} CC) 


Figure 3-6. Comparison of experimental and numerical results of (a) heat flow of concrete 
specimen containing paraffin oil (PCM) in the LGCC experiment, (b) exothermic peak 
representing the heat release of PCM during the solidification process, (c) endothermic peak 
representing the heat absorption of PCM during the melting process. 
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Figure 3-7 (a) Schematic configuration of the slab tests with adapted finite-difference grids, and 
(b) temperature profile at the top of the slab test, T (x=1, t). 
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The nature of heat transfer between the top of the slab and the environment can be classified 
as forced convection as the heat flow is caused by the external means of a fan to change the 
temperature of the environment. The heat transfer coefficient for forced convection using air 1s 
recommended to be within the range of 10 and 250 [W/(m?-K)] (Incropera et al., 2007; Rohsenow 
et al., 1998). Therefore, a heat flux boundary condition is imposed on the top surface of the slab 


(x = 2 mm) to include the effects of convective heat transfer, as described in equation 3-14. 


oT (x, t) 7 


—k, U =h . [r(x = mm, t)- T (x, t)| 3-14 


conv—diss 


Similar to the LGCC test, a difference between the measured heat release during solid PCM 
formation and the associated enthalpy of fusion of paraffin oil due to experimental imperfection 
conditions (thermal bridges, heat leaks, etc.) is observed (Hadi S. Esmaeeli et al., 2017). Hence, a 
sensitivity analysis is conducted for the slab test to obtain the optimal values of input variables; 
while the thermal properties of liquid and solid PCM used in the LGCC and slab tests are assumed 
to be identical, as outlined in Table 3-4. The dissipation coefficients Aconv-aiss and Niatent-loss are 


adjusted to be 20 [WAm7-K)] and 60 %, respectively. 


Table 3-4. Thermal properties of constituents of mortar specimen in the slab test 


Material k [WAm-K)] p (kg/m’) C [KJXke-K] 
Dry concrete 1.75 2250 850 
Liquid paraffin 0.21 670 3071 
Solid paraffin 0.21 800 2461 


The numerical result for temperature profile within a concrete slab without PCM, the 
‘reference’ slab, exposed to a thermal cycle is compared to the experimental data, as shown in 
Figure 3-8. The profile for ambient temperature is shown in Figure 3-8(a) (similar to Figure 3-7(b)) 
to display the cooling and heating moments in the concrete slab. Once the cooling cycle ceases, 
the ambient temperature remains unchanged at -10 °C for 24 h to allow the slab to reach thermal 
equilibrium before the heating cycle begins. Similarly, the ambient temperature remains 
unchanged at +10 °C for 48 h once the heating cycle is completed. 

Figure 3-8(b) and Figure 3-8(c) illustrate the temperature contour within the slab obtained 


from experimental data and numerical results by interpolation of recorded and predicted 
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temperature values at the eight locations within the slab, respectively. In these figures, the white 
color represents the locations within the slab in which the temperature drops below 0 °C. 
Numerical results predict that the temperature within the slab decreases or increases very rapidly 
and is associated with the decrease or increase of ambient temperature in the absence of PCM, 
which is consistent with the experimental results. The absolute temperature difference between the 
predicted numerical results and experimental data is quantitatively calculated at two locations of 
x = 85 mm, and x = 113 mm. It is observed that this difference remains below 1 °C at the location 
of x = 85 mm throughout the thermal cycle. For the location of x = 113 mm, this difference exceeds 
1.2 °C only for a very short period (less than 6 min) and remains below 1 °C for the remainder of 
the thermal cycle. The numerical model also indicates that the temperature within the concrete slab 


without PCM remains below 0 °C for an average time of 25 h. 
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(a) 
Figure 3-8. Thermal analysis of the concrete slab without PCM (a) ambient temperature applied 
to the slab as one thermal cycle, (b) experimental results for temperature contour within the slab 
as a function of time, and (c) numerical results for temperature contour within the slab as a 
function of time. 
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Figure 3-8 continued 
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Figure 3-9 shows the experimental data and numerical results for temperature profile 
within a concrete slab containing PCM exposed to a thermal cycle. The profile for ambient 
temperature is shown in Figure 3-9(a) (similar to Figure 3-7(b)) to display the cooling and heating 
moments in the concrete slab. Figure 3-9(b) and Figure 3-9(c) illustrate the temperature contour 
within the slab obtained from experimental data and numerical results by interpolation of recorded 
and predicted temperature values at the eight locations within the slab, respectively. The analysis 
of experimental data and numerical results shows a relatively slow variation of temperature as the 


gradual heat release during freezing of PCM, or gradual heat absorption during melting of PCM 
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keeps the temperature within slab above 0 °C, for a longer time compared to the temperature within 
reference slab. It is concluded that the gradual heat release or absorption can compensate for the 
heat loss or gain in the slab as the ambient temperature varies. The numerical model implies that 
the slab with PCM remains “warmer” for a longer period; the time at which the temperature is kept 
above 0 °C is reduced to 15 h. Although no ice formation within concrete pores in the slab test was 
considered, the slab would be subjected to freeze-thaw cycling at a temperature below 0 °C, if the 
slab was fully saturated with water solution (1.e., vw = 100 %) and no undercooling would occur. 
Therefore, the incorporation of PCM would reduce the freeze-thaw time by 9 h, and thus improve 
the freeze-thaw performance of the slab. 

Similar to the analysis of reference slab, the absolute temperature difference between the 
predicted numerical results and experimental data is quantitatively calculated at two locations of 
x = 85 mm, and x = 113 mm. It is observed that this difference always remains below 1 °C at these 
two locations throughout the thermal cycle. It is concluded that the numerical model reasonably 
predicts the thermal behavior of concrete slabs with and without PCM exposed to one thermal 


cycle. 
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(a) 
Figure 3-9. Thermal analysis of the concrete slab with PCM (a) ambient temperature applied to 
the slab as one thermal cycle, (b) experimental results for temperature contour within the slab 
as a function of time, and (c) numerical results for temperature contour within the slab as a 
function of time. 


78 


Figure 3-9 continued 
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3.6.4 Field application 


Although the thermal diffusivity (1.e., thermal conductivity divided by heat capacity) of 
PCM incorporated within the concrete infrastructure, such as bridge decks, was determined to be 
an important property governing the effectiveness of PCM at improving the freeze-thaw 
performance (A.R. Sakulich & Bentz, 2012), little is known about the contribution of the PCM 


transition temperature. According to the modeling, the comparison of temperature contours within 
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the concrete slabs with and without PCM indicated that incorporation of PCM with an initial 
freezing temperature above 0 °C significantly reduces the time at which the slab experiences a 
temperature below 0 °C. It is believed that a slight change in the transition temperatures (1.e., the 
rate of heat release or absorption over the transition temperature range) of PCM would 
significantly impact the thermal behavior of the concrete slab exposed to thermal cycles. Therefore, 
three different PCMs are proposed to be used: a paraffin oil, FT4-MT5, and two modified paraffin 
oils, FTO-MT1, and FT2-MT3. FT4-MT5 represents the paraffin oil used in the LGCC and slab 
tests, while initial freezing and melting temperatures of FTO-MT1 and FI2-MT3 PCM are 


modified, as outlined in Table 3-5. 


Table 3-5. Transition temperatures and effectiveness of three PCMs used in field application 


study. 


PCM FTO-MT1 FT2-MT3 FT2-MT3 
Tr (°C) 0 -10 -8 

Tr (°C) 2 -10 -8 

Tmi (°C) 4 -10 -8 

Tu’ (°C) 1.1 3.1 5.1 


Locations where effective more than 10 % 


and less than 20 % (#) = 2 oe 
Locations where effective more than 20 % 37 43 14 
and less than 30 % (#) 

Locations where effective more than 30 % 57 5 2 


Figure 3-10(a) shows the correlation between AH F and freezing temperature for the three 
PCMs used in the numerical model to simulate the gradual heat release due to PCM solidification. 
Figure 3-10(b) shows the correlation between AH F and melting temperature for the three PCMs 


used in the numerical model to simulate the gradual heat absorption due to PCM melting. 
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Figure 3-10. The calculation of modified fusion of enthalpy, AH F as a function of temperature 


during (a) freezing of PCM in the cooling cycle, and (b) melting of PCM in the heating cycle. 


The numerical model is used to predict the thermal behavior of reference pavement, and 
pavements with PCMs incorporated with a dosage of 170 (kg/m’) at 210 different cities throughout 
the United States. The concrete pavements are subjected to arbitrary realistic climate data, 
including hourly averages of temperature for one year from August 2016 through August 2017 
(obtained from the national climate database (Https://Www. Wunderground.Com/, n.d.)) for each 
location. The year-long temperature data is used to simulate the thermal cycling condition imposed 
on the top surface of the pavement to estimate the PCM effectiveness in determining the long-term 
response of the pavement concerning reducing the number of freeze-thaw cycles. For this analysis, 
the only one-year analysis was done. Similar to the thermal conditions imposed in the slab test, 
these pavements are assumed as semi-infinite media in which: (1) one-dimensional heat conduction 
determines the temperature profile through the depth of the pavement, (11) all the thermal properties 
become temperature-dependent when PCM experiences phase transition, (111) the convection 
coefficient h, at the top surface of pavement remains constant, (iv) radiation exchange between the 
pavement and sky is considered to be obstructed by cloud-covered condition, (v) no precipitation 
event is considered to occur. Therefore, the pavement surface conditions used in the model have 
been restricted to those found in the slab test, which includes a dry surface (1.e., free of snow, ice 


and rain), totally cloud-covered, and windy surface. It 1s necessary to note that the equivalent 


81 


temperature of radiation of the sky is logically close to the ambient temperature under cloud- 
covered conditions (Mammeri et al., 2015; Walker & Anderson, 2016). 

These locations are assumed to be fully saturated with water solution; thus, it is subjected 
to freeze-thaw cycling at a temperature below 0 °C. However, the freezing of water solution inside 
the concrete pores occurs at a temperature between ~ -3 to -6 °C (W. Li et al., 2012) due to the 
presence of ions in the pore solution. For the sake of simplicity, this study assumes that the concrete 
pavements experience one freeze-thaw cycle once the pavement temperature reduces to 0 °C. The 
temperature profile within a pavement section is predicted based on the assumed environmental 
inputs on the ambient and the concrete material properties; the numerical model predicts the 
concrete surface temperature (see equation 3-14)) (Antaki, 1997; D. Wang & Roesler, 2012). The 
temperature of the surface of concrete, however, may be different than the air temperature due to 
different factors such as solar radiation, wind velocity, and earth thermal effect that can be the 
objectives of future research. Two criteria are thus defined to determine the PCM effectiveness at 
improving the freeze-thaw performance by comparing the pavement containing PCM with 
reference pavement. These two criteria are: (1) the percentage difference of maximum depth 
within the pavement that experience temperature below 0 °C as shown in equation (3-15a), and (2) 
the percentage difference of time within the pavement that experience temperature below 0 °C as 


shown in equation (3-15b). 


max(X „e [T < OC C),t]) — max(X poy [T < OF C),t]) 
eee ES X 


100 (3-15a) 
max(X „ep [T < O° C),t]) 


X PCM-eff (%)= 


tref [T < OC C), x = 2mm] —t pem [T < OCC), x = 2mm] 


, x100 (3-15b) 
tref [T < 0C C), x = 2mm] 


t pcm -epf (7) = 


where Xpcm-ep1s the percentage difference of maximum depth within the pavement (%), tPcm-eff1s 
the percentage difference of time within the pavement (%), Xref 1s the depth of reference pavement 
which experiences a freeze-thaw cycle at time t, Xpcm is the depth of pavement with PCM which 
experiences a freeze-thaw cycle at time t, tref is the number hours that reference pavement 
experiences a freeze-thaw cycle on the top surface (x = 2 mm), tpcu 1s the number hours that 


pavement with PCM experiences a freeze-thaw cycle on the top surface (x = 2 mm). 
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Modeling of reference pavement shows that there is a total of 44 locations where PCM 
effectiveness is less than 10 %, as these locations are either warm enough that freezing cycles 
rarely occur (located primarily in the Southwest, Southeast and along the Gulf Coast), or cold 
enough that the temperature of the pavement generally remains below 0 °C (located primarily in 
the North and the Northeast). Figure 3-11 and Table 3-5 show the effectiveness of the three PCMs 
on the freezing time and depth of concrete pavements throughout the United States. It is observed 
that all three PCMs are effective in regions with milder weather (1.e., the majority of cities located 
in the East, Central, and Northwest regions). In a total of 78 cities, one of the three PCMs reduces 
both freezing time and depth of pavement by at least 10 percent. Sakulich et al. also predicted that 
the incorporation of two PCMs with transition temperatures of 4 °C and 6°C within concrete 
mixtures was effective in these regions, where at least 48 of the 100 most populous American cities 


are located (A.R. Sakulich & Bentz, 2012). 
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(a) 
Figure 3-11. Map of cities in which the concrete pavements incorporating three PCMs investigated 
here decrease the time and depth of freezing. The types of PCM used in pavements are (a) FTO- 
MT 1, (b) FT2-MT3, and (c) FT4-MTS5. Point color indicates what percentage the time and depth 
of freezing within concrete pavement decreased, by quartile: blue indicates PCM effectiveness 
less 10 %, green indicates PCM effectiveness between 10 % and 20 %, orange indicates PCM 
effectiveness between 20 % and 30 %, and red indicates PCM effectiveness more than 30 %. 
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Figure 3-11 continued 
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In this study, the numerical model predicts that PCM incorporation is more effective at preventing 
freeze-thaw cycling, as the freezing and melting temperatures are decreased toward 0 °C for FTO- 


MTI and FT2-MT3 PCMs. Figure 3-11(a) shows that the FTO-MT1I PCM is relatively more 
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effective in reducing the freezing time and depth in states with milder weather compared to the 


effectiveness of FT2-MT3 and FT4-MT5S PCMs, as displayed in Figure 3-11(b) and Figure 3-11(c), 


respectively. It is concluded that the incorporation of a PCM with relatively low ie would lead to 


release a considerable amount of heat inside the pavement right before the ice formation within 


concrete pores and prevent localization of freeze-thaw damage. Conversely, a lower T% would 


cause the solid PCM to melt completely to retrieve the heat right before a subsequent cooling cycle, 


thereby decreasing the localized damage due to freeze-thaw cycles. 


3.7 Conclusions and summary 


This paper investigated the development of a computer simulation tool that used a one-dimensional 
finite difference method and homogenization techniques to assess the thermal behavior of 
cementitious materials containing PCM exposed to thermal cycling. The computational results 
were compared to the experimental data for mortar specimen tested in the LGCC (small-scale) test 
with a temperature range of -40°C to +24°C and concrete slabs (large-scale) tested with a 
temperature range of -10 °C to +10 °C. The numerical model was validated by measured data from 
LGCC experimental data where it was concluded that the consideration of a gradual phase 
transition of PCM with a narrow range of transition temperature better predicted the thermal 
response (temperature profile and associated heat flow) of the specimen, while the LWA pore sizes 
had a negligible effect. Analysis of the slab test also showed that the time of freezing at which 
temperature reduced below 0 °C was decreased roughly by 9 h when compared to a slab without 
PCM, which implied the improvement in the freeze-thaw performance of slab. 

The numerical model was also used to investigate the PCM effectiveness on reducing the 
impact of freeze-thaw cycling within concrete pavements located in different regions of the United 
States. It was observed that the transition temperatures of PCM during freezing and melting events 
are important properties that govern the effectiveness of PCM at reducing the impact of freeze- 
thaw cycles. In other words, they help control how quickly the temperature of the pavement drops 
below the freezing point of pore solution within concrete pores. In 122 out of 210 cities 
investigated, a dose of 170 (kg/m’) of FTO-MT1 PCM would reduce the freezing time and depth 
of concrete pavement by at least 10 %, primarily throughout the East, the Central, and Pacific 


Northwest regions. The model results also showed that low transition temperatures were more 
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effective, where pavements located in an additional 13 and 44 experienced a reduction in freezing 
time and depth of pavements by at least 10 %. These preliminary results suggested that the use of 
PCM-LWA composites as a thermal energy storage system shows promise in improving the 


freeze-thaw performance of concrete pavements and merits further, more detailed investigations. 
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4. A TWO-STEP MULTISCALE MODEL TO PREDICT EARLY AGE 
STRENGTH DEVELOPMENT OF CEMENTITIOUS COMPOSITES 
CONSIDERING COMPETING FRACTURE MECHANISMS 


This chapter contains work that was originally published in Elsevier as “Hadi S. Esmaeeli, 
Mehdi Shishehbor, W. Jason Weiss, and Pablo D. Zavattier1. A two-step multiscale model to 
predict early age strength development of cementitious composites considering competing fracture 
mechanisms. In Journal of Construction and Building Materials. Volume 208., pp. 577-600. 
Elsevier, 2019.” The original article has been used with permission as stated below. 

Reprinted by permission from Elsevier, in: Journal of Construction and Building Materials. 
“A two-step multiscale model to predict early age strength development of cementitious 
composites considering competing fracture mechanisms.” Hadi S. Esmaeeli, Mehdi Shishehbor, 
W. Jason Weiss, and Pablo D. Zavattieri, COPYRIGHT (2019). 


4.1 Introduction 


Mesomechanical models have recently gained prominence in the analysis of damage and 
failure mechanism of concrete (Benedetto et al., 2018; Dutta & Chandra Kishen, 2018; Elias et al., 
2015; Zhichao Liu & Hansen, 2016; Lopez et al., 2008b; R. Zhou & Lu, 2018). A well-recognized 
approach to fracture for concrete is to use the cohesive zone model which enables the crack to 
propagate through zero-thickness interfaces equipped with a normal-shear traction-separation law 
representing non-linear fracture process (Feng & Wu, 2018; Toyama et al., 2018; X. Xu & 
Needleman, 1995). Several FEM-based models in the literature have dealt with cracking in 
heterogeneous materials via computational models, such as continuum damage mechanics (CDM) 
(Kim & Abu Al-Rub, 2011), extended finite element method (XFEM) (Ng & Dai, 2011) and 
cohesive zone model (CZM) (W. Gao et al., 2015). CDM model cannot model the initiation, and 
propagation of explicitly, as cracks occur continuously. In the XFEM approach, additional degrees 
of freedom are considered in the standard finite element using local enrichment functions, where 
further external criteria are required to anticipate the cracking (Sukumar et al., 2000). A major 
drawback of this method is the complexity of problems involving pervasive fracture and 
fragmentation (1.e., interaction of many cracks which 1s of interest in this paper). On the other hand, 
the cohesive element approach is consisted in inserting interface (cohesive) elements along the 
surfaces of the two- or three-dimensional continuum elements correspondingly (Camacho & Ortiz, 
1996; M. Ortiz & Pandolfi, 1999; Rimoli & Rojas, 2015; X.-P. Xu & Needleman, 1996). Finite 


element method (FEM) based models have been beneficial to analyze the mesostructure and the 
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effects of local properties on the macroscopic behavior of concrete. For instance, Wang et al. (X. 
Wang et al., 2015) and Lopez et al. (Lopez et al., 2008a) employed a mesomechanical model to 
systematically study by considering constant cohesive properties (strength and fracture energy per 
unit area) of mortar and ITZ* phases and morphological properties (arrangement, size, and volume 
fraction) of coarse aggregate phases. The mechanical response was calculated to consider 
microcracking coalescence, localization, and macroscopic crack formation. While they did not 
consider the aging effect for the cement paste and local failure of aggregate phase on early-age 
mechanical response of concrete, their model provided information on how cohesive properties of 
concrete constituents have direct effect on the resulting strength, softening behavior, and crack 
morphology. 

This paper proposes a computational approach to assess the competing mechanisms 
occurring under remote predominantly mode-I failure for early-age mortar and concrete. The study 
analyzes the evolution of the mechanical properties for the cement paste and ITZ phases, along 
with constant mechanical properties for the aggregate phase to simulate the crack propagation 
through the potential phases to complete fracture. We hypothesize that the knee point is primarily 
controlled by the competition between the fracture of the cement that change over time at early 
ages and the fracture of the aggregate phase at later ages. This work aims to develop a holistic 
numerical framework to understand the role played by the aggregate phase on limiting the early- 
age mechanical response and resulting alteration in the cracking pattern of cementitious 
composites. 

The multiscale approach consists of a two-step homogenization model where (1) the material 
properties of the cement paste, fine aggregate, and ITZ™ (“m” denotes mortar) phases are utilized 
to analyze the homogenized properties and associated failure mechanisms for heterogeneous 
mortar using a coupled continuous and discontinuous homogenization scheme (De Borst, 2008; 
Geers et al., 2010; Nguyen, Stroeven, et al., 2011; Phu Nguyen et al., 2010; Verhoosel et al., 2010), 
and compares to the experimental observations, and afterward (2) the homogenized mechanical 
properties of mortar obtained from the first step are upscaled to be utilized along with the 
mechanical and morphological properties of coarse aggregates to quantify the mechanical response 
and crack propagation within mesostructure of concrete via continuous homogenization scheme 


and contrast them with experimental counterparts. 
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The proposed multiscale approach considers three representative structures of cementitious 
composites, outlined in Figure 4-1 (Bernard et al., 2003; Constantinides & Ulm, 2004): 
Level I (10-104 m): The microstructure of paste material with the characteristic length of lp is 
composed of large calcium hydroxide (CH) crystals (<200 um), reacted aluminate phase (A-phase), 
pore space, and calcium silicate hydrate (C-S-H) (Bernard et al., 2003; Constantinides & Ulm, 
2004). The structure of paste develops during the hydration process, which leads to reducing the 
porosity and increasing the strength (J. J. Thomas & Jennings, 2006). The size of homogeneous 
phases is the order of magnitudes smaller than the sizes of mortar and concrete specimens used in 
four-point flexural tests (Barde et al., 2005); thus, the paste material is treated as a homogeneous 
and brittle phase in this study. 

Level II (10°-107 m): The mortar material is treated as a three-phase composite with the 
characteristic length of /m that is composed of homogeneous phases of paste as the matrix, and fine 
aggregate particles as inclusions (<4.75 mm), and ITZ™ embedded between cement paste and fine 
aggregate. This level has been the focus of previous mesomechanical models, both numerically 
(G. Li et al., 1999) and analytically (Garboczi, 1993). In this work, Level II is considered as the 
starting point for the computational homogenization approach. 

Level IH (10-10! m): The concrete material is considered as a three-phase composite with 
the characteristic length of le that is composed of homogeneous phases of mortar as the matrix, 
coarse aggregate as the inclusion (4.75 mm< <25 mm), and ITZ* (“c” denotes concrete) embedded 


between the mortar and coarse aggregate particles (Garboczi, 1993; G. Li et al., 1999). 
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Figure 4-1. The three-level representative structures of cementitious composites, adapted from 
(Bernard et al., 2003; Constantinides & Ulm, 2004). 


4.2 A two-step homogenization mesomechanical model 


A mesomechanical model has been developed to understand better the composite effect on 
the early-age mechanical response of cementitious composites under mode-I loading. This 
includes the prediction of tensile stiffness and strength gain, as well as damage mechanism through 
micro-crack propagation and coalescence, macro-crack pattern, and fracture energy. Finite element 
method (FEM) based models have been beneficial to analyze the mesostructure and the effects of 
local properties on the macroscopic behavior of concrete. For instance, Wang et al. (X. Wang et 
al., 2015) and Lopez et al. (López et al., 2008a) employed a mesomechanical model to 
systematically study by considering constant cohesive properties (strength and fracture energy per 
unit area) of mortar and ITZ* phases and morphological properties (arrangement, size, and volume 
fraction) of coarse aggregate phases. The mechanical response was calculated to consider 


microcracking coalescence, localization, and macroscopic crack formation. While they did not 


90 


consider the aging effect for the cement paste and local failure of the aggregate phase on early-age 
mechanical response of concrete, their model provided information on how cohesive properties of 
concrete constituents have a direct effect on the resulting strength, softening behavior, and crack 
morphology. 

This work develops a two-dimensional FEM-based mesomechanical model that accounts for 
the evolution of tensile stiffness and strength of cement paste (and ITZ) in time, and the differences 
between the strengths of constituents. This study also investigates the effects of Poisson’s ratio 
and fracture energy development of cement paste (and ITZ) on elastic and fracture properties of 
cementitious composites. To test the validity of the hypothesis developed in the previous section, 
the mesostructures of mortar and concrete are modeled under uniaxial tension and four-point 
flexural tests. Prior experimental observations illustrate that different failure mechanisms are at 
work as the hydration of cement paste progresses. We hypothesize that the effects of cement paste 
strength versus aggregate strength compete in the hydration process: before knee point, the tensile 
strength of cement paste phase seems to have a stronger influence than the tensile strength of 
aggregate phase on the evolution of the composite strength of cementitious composites, and after 
knee point, it is the inverse. As the cement paste gains strength, it also becomes more brittle, which 
forces the cement paste to fail primarily in tension. The failure mechanism is observed to be 
dictated by crack propagation around the aggregates at a low degree of hydration. At the same 
time, the knee point in the bilinear mortar and concrete curves represents a transition to aggregate 
fracture. 

To provide insights into the mechanical properties of mesostructure constituents and how 
different phases fail during mode-I loading, we explicitly consider the competition of the crack 
propagation of each phase in our mesomechanical model. The contribution of the aggregate phase 
of different sizes to the macroscopic mechanical response and cracking patterns occurs at two 
length scales. The two levels that particularly distinguish the contribution of fine aggregates at 
Level II from the contribution of coarse aggregates at Level III indicates a two-step 
homogenization scheme, as shown in Figure 4-2. The homogenization procedure is used as a 
computational tool to alleviate the burden of material properties characterization, in which the 
properties can be homogenized as a function of key morphological and mechanical properties such 


as the structural composition, aggregate morphology, and the mechanical properties of constituents. 
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Figure 4-2(a-c) illustrates a continuum-based finite element modeling framework that 
considers the individual mechanical properties of the mortar constituents (e.g., cement paste and 
fine aggregate) to homogenize the elastic and cohesive properties of Level II mortar to replace the 
heterogeneous mortar by an identical homogeneous mortar in Level III concrete. To this end, we 
employ continuous and discontinuous computational homogenization schemes. Figure 4-2(a) 
shows an arbitrary mesostructure of heterogeneous mortar. To characterize the linear elastic region 
of mortar, the developing tensile stiffness of cement paste E” (t) over early ages were obtained from 
experimental measurements (Yoshitake et al., 2012; Zhao et al., 2014) along with constant tensile 
stiffness of fine aggregate (BE), as displayed in Figure 4-2(b). The inset in Figure 4-2(b) shows the 


evolution of Æ” as a function of t (e.g., E- > E- where t > t), and a constant E which is 


determined experimentally by Esmaeeli et al. (Hadi Shagerdi Esmaeeli, 2015) based on 
mechanical analysis of limestone beams under three-point flexural test. Afterward, the tensile 
strength development of mortar f/"(t) is determined by specifying the cohesive properties (tensile 


strength and fracture energy) of the interface elements embedded within cement paste (f/(t) and 


GÈ), fine aggregates (f/“ and Go and ITZ™ (f(t) and GÈ) continuum elements, as illustrated in 
Figure 4-2(c). We formulate the traction-separation law (T-u curve) for early-age mortar at each 
discrete time (f/”(f) and Gi) to capture the nonlinear pre-peak region governed by tensile micro- 
cracking and post-peak region dominated by micro-crack coalescence, macro-crack formation, and 
corresponding shear failure development. The interpretation of this T-u curve is demonstrated in 
the inset in Figure 4-2(c). T-u curve for cement paste evolves as its tensile strength develops 
(e.g., A > fe where t2>t7) acquired from experimental data (Barde et al., 2005), while GÈ. 1S 
considered to remain constant. For fine aggregate, T-u curve is characterized by constant f/“ and 
GÈ. obtained from experimental results (Hadi Shagerdi Esmaeeli, 2015). We will get back to this 
nonlinear fracture approach associated with a bilinear traction-separation law. 

The second homogenization step is composed of connecting from Level II to Level HI, where 
the continuous homogenization scheme is used to quantify the macroscopic tensile stiffness and 
strength of early-age concrete (f(t) and E‘(t)). Figure 4-2(d) shows an arbitrary mesostructure of 
heterogeneous concrete composed of the mortar matrix resulting from the first homogenization 
step and coarse aggregate particles. Similar to the modelling of Level II mortar, we adopt a 


continuum FEM-based model to characterize the linear elastic region of Level III concrete by 
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utilization of developing tensile stiffness of mortar at discrete ages (Ez) > Ey) where t2 > t7) along 
with constant tensile stiffness of coarse aggregate (E), displayed in Figure 4-2(e). Unlike Level 
II, here we employ the values for properties E(t), f/"(t), and G; taken from first homogenization 
step, while the properties E“, ff“, and Gf are obtained from experimental measurements. Similar 


to Level II, Figure 4-2(f) depicts the evolution of T-u curves for mortar at early ages as fét, > fre, 


where f2>t;, while Gi remains constant. The mixture proportioning and specimen preparation 
procedures for making the mortar and concrete specimens are explained in Appendix A. It is 
necessary to note that the viscoelastic phenomenon of cement paste is neglected in this study, as 
the mortar and concrete specimens are tested at fast loading rate. A few mesoscopic models have 
systematically investigated the effects of temperature and moisture variations on the mechanism 
of shrinkage cracking in cementitious composites (Neubauer et al., 1996; Tang et al., 2013, 2016). 
This study, however, is focused on well cured and unrestrained mortar and concrete specimens, 
where the effect of shrinkage attributed to moisture diffusion due to external drying does not 
contribute to the cracking. Future work will focus on better determining and quantifying if 


moisture diffusion results in the development of stress and early-age cracking. 
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Figure 4-2. Numerical framework representing a two-step multiscale modeling the mechanical 
response development of (a) a heterogeneous Level IJ mortar, knowing (b) the tensile stiffness of 
paste and fine aggregate (E’(t) and E), and (c) the tensile strength (fP (A) and ff“) and fracture 
energy (Gand GÍ”) of paste and fine aggregate; afterwards modeling the mechanical response 
development of (d) a heterogeneous Level III concrete, knowing (e) the tensile stiffness of mortar 
and coarse aggregate (E”(t) and E“), and (P) the tensile strength (f/"(t) and f£*) and fracture energy 
(Gj and G°) of mortar and coarse aggregate. 
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After finite element meshing of mesostructures of cementitious composites”, 4-node zero- 
thickness cohesive interface elements are pre-inserted between existing 3-node bulk elements by 
developing an in-house code in C programming language. We employed an explicit finite element 
approach enhanced with a central difference time integration scheme to predict the behavior of 
heterogeneous mesostructures in transient dynamics applications using finite element package 
FEAP (Taylor, 2008). To calculate the field variables at respective nodes, the requirement for 


small time steps for stability reason is met. 


4.2.1 Mesomechanical representation 
4.2.1.1 Representative volume element (RVE) 


Homogenization methods have been used to estimate the macroscopic material properties 
by averaging the response of a representative volume element (RVE) with a determined size where 
the mesostructure 1s specifically analyzed (Nguyen, Stroeven, et al., 2011). For materials having a 
random mesostructure, this RVE has been used for the description of composition and 
heterogeneity of the internal structure of the material system (see Figure 4-2(a)). The RVE 
describes the different phases of mesostructure with their specific structural geometry, in which 
the spatial distribution of material properties has been addressed. Considering the development of 
mechanical properties of matrix and ITZ phases at an early age, we incorporate these time- 
dependent properties into separate RVEs to perform the simulations at discrete ages. It is worth 
mentioning that an RVE represents a mesoscopic specimen statistically when (1) the homogenized 
properties do not alter considerably as the size of specimen increases, and (2) the homogenized 
properties do not depend on the mesostructural randomness as the size of the specimen is large 


enough (Nguyen, Lloberas-Valls, et al., 2011). 


4.2.1.2 Fracture-based interface cohesive law 


A mesomechanical model has been developed to assess the roles of elastic deformation of 
continuum elements coupled with the damage and failure of the material on the mechanical 
behavior of heterogeneous RVEs selected for cementitious composites. The existence of phases 


with dissimilar properties requires the establishment of an approach to fracture, where the 


? Commercial FE package ABAQUS is used to generate triangular bulk elements with random sizes and orientations. 
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contribution of geometrical and mechanical properties of heterogeneities on the overall mechanical 
behavior of cementitious composites is determined. While previous FEM-based models were 
beneficial to quantify the stress distribution and local material failure in the mortar and interface 
between mortar and coarse aggregates (ITZc) (Mirzabozorg & Ghaemian, 2005; Skrikerud & 
Bachmann, 1986), not many have addressed the influence of mechanical properties of hydrating 
cement paste and aggregates resulting in potential cracking in the aggregates. By inserting 
interface elements between continuum elements, these elements represent the potential crack paths 
that allow the cracks to propagate through them. Several FEM-based models in the literature have 
dealt with cracking in heterogeneous materials via computational models, such as continuum 
damage mechanics (CDM) (Kim & Abu Al-Rub, 2011), extended finite element method (XFEM) 
(Ng & Dai, 2011), and cohesive zone model (CZM) (W. Gao et al., 2015). CDM model cannot 
model the initiation, and propagation of explicitly, as cracks occur continuously. In the XFEM 
approach, additional degrees of freedom are considered in the standard finite element using local 
enrichment functions, where further external criteria are required to anticipate the cracking 
(Sukumar et al., 2000). A major drawback of this method is the complexity of problems involving 
pervasive fracture and fragmentation (1.e., interaction of many cracks which is of interest in this 
paper). On the other hand, the cohesive element approach is consisted of inserting interface 
(cohesive) elements along the surfaces of the two- or three-dimensional continuum elements 
correspondingly. This approach suits for the analysis of problems which include pre-defined crack 
paths; however, several recognized issues have effects on its accuracy when modeling the 
problems with random crack directions, 1.e., (1) artificial flexibility resulting in alteration of 
propagation of stress wave speed, (11) spurious crack tip speed, and (111) mesh size and orientation 
effects. Regardless of these widely definite constraints, this approach is one of the most robust 
methods for fracture analysis. Therefore, we employ a non-linear fracture analysis in which the 
zero-thickness interface elements enhanced with a traction-separation cohesive law embedded as 
intra-phase (aggregate-aggregate and matrix-matrix) and inter-phase (aggregate-matrix) elements. 

It is assumed that the interface elements carry forces that oppose separation and shear 
between two surfaces until failure, where the magnitude of these forces is a function of the relative 
separation and shear displacement between the two surfaces. Note that the intrinsic strength and 
fracture toughness of cementitious composites are much larger in pure mode-II than those in pure 


mode-I. Therefore, the tensile stresses at the intra- and inter-phase elements have a significant 
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influence on the load capacity and fracture behavior of these materials subjected to remote mode- 
I loading condition. In contrast, the contribution of shear tractions at these interface elements to 
the overall mechanical behavior is negligible. The formulation of a cohesive element approach as 
a function of traction (rather than strain) and Ge regularizes unstable behavior during crack 
nucleation occurring on the weak interface element. The continuum elements are commonly 
considered to show a linear elastic behavior. For instance, let us consider a 2D Level II mortar 
RVE discretized with linear triangular elements, shown in Figure 4-3(a), with prescribed 
displacement, u on the Dirichlet boundary condition ’. Schematics of three sets of interface 
elements with different cohesive laws inserted within fine aggregate particles (fa-fa), cement paste 
matrix (p-p), and their interface (ITZ™) is shown in Figure 4-3(b). While characterizing the fracture 
properties of ITZ phase at two levels of H (TZ™ in the mortar) and III (ITZ* in concrete) can be 
useful to understand the onset of cracking as a displacement discontinuity between matrix and 
particles with dissimilar properties may exist, we take an approach in this study where the cohesive 


properties of both ITZ phases are explicitly considered to be identical to their associated matrix. 


Continuum 
Elements 


E p 
E ta 


Interface 
Elements 


a p-p 
|] fa - fa 
E ITZ™ 





(a) (b) 


Figure 4-3. (a) Fracture analysis of an RVE for Level IJ mortar RVE subjected to mode-I loading 
containing (b) continuum elements for two phases of the matrix, and aggregate as well as three 
interface elements inserted between paste-paste (p-p), fine aggregate-fine aggregate (fa-fa), and 
paste-fine aggregate (ITZ™), with a cohesive law defined in (c) opening mode, and (d) shear 
mode. 
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Figure 4-3 continued 
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The cohesive law for interface elements is formulated in terms of normal and shear 
components of stresses o = (Tan, T;) on the interface element and corresponding relative 
displacements u = (Un, ut). We employ a bi-linear cohesive law in opening and shear mode 


implemented as a user-subroutine in the finite element package FEAP. The non-dimensional 


2 2 
displacement jump, À = (=) + (=) is defined to describe the degradation of interface 
n t 


element by normalizing the evolution of un and ur parameters respect to the maximum normal and 
shear displacement values to reach failure (ôn, 0+). 

The mechanical behavior of the interface elements is primarily governed by two fracture 
energy parameters Gr (associated with opening mode) and Grr (associated with shear mode). Two 
traction components describe the resistance to separate the interface elements in normal (T, (t) = 


1—A* =n) Tmax (t) : Ee (=) Tmax(t) 
r ( 5,) Ane ) and tangential directions (7;(t) = es) ie 





), Where Tmax and Tmax are 





nominal tensile and shear strength. The value of 2° monotonically increases given by 4*= max (4*, 
A), where Amax = Acr initially and Amax = 2 1f Amax > 2. It should be noted that the normal and tangential 
tractions reach their maximum value at A = Acr. The cohesive law in opening and shear modes is 
identified in Figure 4-3(c) and (d), respectively. It should be noted that the initial stiffness of 


cohesive law is selected with care to avoid the addition of compliance due to the presence of 
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interface elements if compared with the compliance of continuum elements. Therefore, at this 
initial stage, the continuum elements remain intact. 

The area under the traction-separation curves corresponds to fracture energy in the opening 
(Gpc=0.5Tinax.0n) and shear modes (Gic=0.5tmax.0r) required to create new discontinuous surfaces 
(i.e., cracks). When displacement initially increases in the linear-elastic portion of the traction— 
separation response, the interface element inserted around the crack tip quickly raises to its strength 
and progresses to the softening portion of the response. Once 1 >/cr, the interface element initiates 
to degrade, thus dissipating energy. Further increasing the tip displacement, elements adjacent to 
the crack tip also undergo reversible deformation, in which the characteristic length of these 
elements is defined as cohesive zone length, lez ~ EG,-/T,Agx. Any unloading during degradation 
follows the irreversible path illustrated in Figure 4-3(c) and (d). The onset of cracking in the 


interface element initiates while Amax=1 (1.e., T,=T;=0). 


4.2.2 A coupled multiscale approach for cohesive crack homogenization 


Figure 4-5 presents the fundamental concept of continuous and discontinuous 
homogenization schemes for cohesive crack modeling. For quasi-brittle materials such as 
cementitious composites subjected to a simple uniaxial tension, the averaged stress-strain response 
(o>€) can be separated into two regions of pre-peak hardening and post-peak softening, as shown 
in Figure 4-5(a). The early part of the hardening region of o—é curve (highlighted with a lighter 
color in Figure 4-5(a) determines the elastic properties (£) and Poisson’s ratio (v) utilizing 
continuous homogenization scheme. The continuous homogenization calculates a homogenized 
macroscopic o—é€ response governing the bulk behavior before strain localization (Geers et al., 
2010; Smit et al., 1998). For completeness, this computational homogenization scheme is 
described here. 

The behavior of bulk material at a macro level is determined through the continuous 
computational homogenization as described in (Zohdi & Wriggers, 2008). For completeness, this 
homogenization scheme is briefly described as follows. The homogenized stresses and strains of 


an RVE are determined as the volume averages of its element constituents using equation 4-1. 


1 1 
o = ales Oe dN = me tj x; ar, 4-] 
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E€ = ! | Ee dQ = | (ufn? +ufnf) dr 
lo, PARA 


where denotes the RVE domain, |.Q,| is the area of RVE, T, is the boundary of RVE, ce and ee 





are the stress and strain tensors at bulk Gauss points, respectively, t represents the traction vector 
while x* is the position vector and n‘ is the unit outward normal to 1}. Therefore, the homogenized 
tensile stiffness, Dois determined using the analytical expression o = (Do.e). 

For modeling cohesive cracks in random heterogeneous quasi-brittle materials, a 
discontinuous homogenization scheme is employed as a regularization technique for the 
continuous homogenization to extract the localized deformation from the state of deformation of 
RVE though failure zone averaging technique. Therefore, we define the homogenized stresses (Tn) 
and strains (€dam) as the volume averages of the stress and strains of RVE elements within 24 rather 
than over © using equation 4-2. 

T, = FA Jy, %e Wes Edam = a fo, Ee Aa 4-2 

The above integrals are calculated directly using numerical quadrature as they cannot be 
converted to surface integrals. We, thus, define un as the damage opening tigg,, shifted to the left 
by Udam to obtain the homogenized traction-separation relations using equation 4-3. 

Un = Udam — Udam, Udam = Edam * (lpn) 4-3 
where n denotes the cohesive crack normal vector. 

To determine the existence of an RVE for softening mortar when using continuous and 
discontinuous homogenization scheme, let us consider the example shown in Figure 4-4. Table 
4-1 outlines the material properties for continuum elements as well as interface elements at two 
early ages of 16 h and 96 h. The equilibrium path is determined through applying a uniaxial 
displacement control with a fixed increment of 2x10° mm. A plane strain condition is assumed. 
The behavior of bulk mortar can be calculated a priori using the continuous homogenization 
scheme. Given the application of the prescribed displacement field on the RVE, the mesoscopic 
boundary value problem (BVP) is computed to quantify the macroscopic stress through the 
calculation of the mesoscopic stress field in the pre-peak region. Furthermore, the discontinuous 
homogenization scheme is considered as an effective complementary method to the continuous 
homogenization scheme to homogenize the cohesive law in random heterogeneous RVEs, thus 
define a unique RVE for softening materials. To validate the objectiveness of a unique RVE for 


mortar w.r.t its size, the response of five mortar RVEs with /m of 6 mm to 15 mm subjected to 
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uniaxial tension is homogenized through discontinuous homogenization scheme, displayed in 
Figure 4-4. It is necessary to note that, although the lez parameter of interface elements inserted 
along aggregate-aggregate continuum elements is very large respect to the RVE size, it is 
suggested to consider the optimum size of the RVE about three times larger than the fine aggregate 
size (Z. Bazant & Oh, 1983). The convergence of Ta-un relation upon increasing the size of RVE 


verifies the existence of an RVE. 


Table 4-1. Material parameters of continuum and interface elements within mortar RVE. 


Mechanical t= 16h t= 96h 
Unit fa 
Property p ITZ™ p ITZ™ 
E [GPa] 6.162 - 16.342 - 31.074 
D [-] 0.2 - 0.2 - 0.15 
Gic [J/m?] 4] 4] 41 41 43.4 
Gic [J/m"] 410 410 410 410 434 
Tmax [MPa] 1.68 1.68 5.862 5.862 3.8 
Tmax [MPa] 2.33 2:93 8.8 8.8 5.7 
Ôn [mm] 0.0488 0.0488 0.013 0.013 0.0013 
Ôt [mm] 0.349 0.349 0.0931 0.0931 0.015 
lez [mm] 89.51 89.51 19.49 19.49 38.103 
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La = 6 mm la = 8 mm bin = 10 mm 


(c) 








la = 12 mm la = 15 mm 


(d) (e) 


Figure 4-4. Mortar RVEs for the case of 55% fine aggregate by volume with sizes of (a) /m = 6mm, 
(b) lIn = 8mm, (c) Im = 10mm, (d) lIn = 12mm, and (e) In = 15mm. 


When strain localizes within the RVE, the homogenized o—é response exhibits a local 
strain-softening behavior. Subsequently, the response obtained with this homogenization approach 
depends on the size of the RVE, the homogenized response becomes more brittle for larger sizes 
of the RVE (e.g. lm, > lm,), shown in Figure 4-5(a) (Gitman et al., 2007). 

While a unique RVE cannot be found using continuous homogenization scheme, as the 


material does not maintain its statistical representativeness at the onset of strain localization, 
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Nguyen et al. (Phu Nguyen et al., 2010) introduced a novel approach utilizing the failure zone 
averaging technique to develop a discontinuous homogenization scheme, in which the cohesive 
law has been homogenized. Using the failure zone averaging approach, the discontinuous 
homogenization scheme explicitly considers the averaging over a localization domain, denoted as 
2u. Therefore, a cohesive law is upscaled by extracting only the localized deformation using the 
discontinuous homogenization scheme. There is no necessity to calculate the bulk stresses outside 
the localization domain as the material is considered to undergo a linear elastic behavior. It should 
be noted that this model assumes the nonlinear part of the pre-peak response (indicated by the 
darker region) to be insignificant. 

The total deformation for RVE A can be described using the following homogenization 
relation: 


A= lemy] Do Tet igam oe 4-4 


where A is the displacement over the entire domain of RVE, the first RHS term demonstrates the 
linear displacement over the domain (2, Udamrepresents the deformation quantified at the onset of 
crack localization in 24, and un is the separation at one cohesive integration point, as shown in 
Figure 4-5(b). In equation 4-4, l» is the width of the localization band, the matrix Do is the 
homogenized tensile stiffness quantified over the domain ©, and Tris the corresponding traction 
at given un. Computation of Ly and Udam 1s performed using the failure zone averaging approach. 
The reader is suggested to refer to the work of Nguyen et al. (Nguyen et al., 2012b) for further 
detailed information on the formulation and computer implementation. 

By comparison of responses obtained from discontinuous homogenization scheme for 
various RVE sizes, we observe a small discrepancy between the peak stress and post-peak 
responses which are associated with the deviation in the hardening response (1.e., the darker region 
in Figure 4-5(a)) which are not considered in the discontinuous homogenization scheme (i.e., the 
pre-peak region is considered to exhibit a linear behavior). Herein, we use five mesostructures with 
randomly distributed fine aggregates for a certain size of RVE to compensate for this discrepancy 
by quantifying the mean value and standard deviation of Tmax. Nguyen et al. (Nguyen et al., 2012b) 
also demonstrated that both continuous and discontinuous homogenization schemes yield the same 
value of peak stress; Therefore, we consider that Tmax can also be obtained using continuous 


homogenization scheme. 


103 








Figure 4-5. Using continuous and discontinuous homogenization schemes to couple macro-and 
meso level models. (a) Altering homogenized o-é relations obtained with continuous 
homogenization scheme, and (b) unique homogenized Tn-un relations obtained with discontinuous 
homogenization scheme. 


The principle of separation of scales is defined by assuming the RVE size at the meso level be 
significantly smaller than the macro-scale over which the loading 1s applied (Nguyen et al., 2012a). 
As aresult, the stress and strain fields at the macro level are uniform over the RVE. It is necessary 
to mention that the minimum distance between coarse aggregates in the RVE of Level III concrete 
is most likely to be small compared to the size of the RVE of Level II mortar. In this work, this 
issue has been resolved as the numerical solutions will demonstrate a good quantitative agreement 
in contrast with experimental observations. It is also worth mentioning that, in the case of material 
discontinuity during the crack initiation, the size of the interface elements should be small 


compared to the lez parameter to obtain a homogenized cohesive law. 
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4.3 Results and discussion 
4.3.1 Step 1: Mortar composite at Level II 


The mortar material at the meso level is considered as a three-phase composite, in which 
the fine aggregate phase with a relatively high stiffness plays the role of inclusion surrounded by 
a weaker ITZ™ embedded in a paste matrix with relatively low stiffness. From a geometrical 
standpoint, it is indicated that the precise geometry of the inclusions is less likely to considerably 
affect the homogenized elastic properties when the inclusions have a relatively higher stiffness 
than the matrix (Simeonov & Ahmad, 1995). However, the exact geometry of the inclusions is of 
importance to quantify the mechanical and cracking responses. Therefore, 2D mortar RVEs 
digitized from realistic images of mortar are generated using computed tomography (Lavergne et 
al., 2015). The first step is composed of the application of continuous and discontinuous 
homogenization schemes to Level II mortar RVE with a characteristic length of Im = 6 — 15 mm. 
The mechanical properties account for the effect of porosity (capillary and air-entrained pores of 
sizes smaller than 100 um), of hydrated products (< 1 um), and of unhydrated cement particles (< 
100 um) present at smaller scales. The RVE is considered to be composed of three phases: (1) 
cement paste of volume Vp, (11) fine aggregate of volume Vy, and (111) the ITZ of zero thickness 


and volume (VITZm = 0). The volume fractions of these phases are described in Equation 4-5. 


fra = VealVni fo = Vp/Vn = 1- fra 4-5 


where fp and ff parameters represent the volume fractions of cement paste and fine aggregate 
phases. It is worth mentioning that the average width of ITZ is of the order of 15 um around each 
aggregate particle (Scrivener et al., 2004), which is 1/7th of average element size used in mortar 
RVE. Considering a constant non-zero width for ITZ may lead to computationally expensive 
models; thus, we assume an ITZ phase with zero width. 

The early-age elastic properties of the three phases in the mortar and the cohesive properties of the 


interface elements are defined as follows: 


e Regarding the Level I paste, the development of tensile stiffness, EP (t) is obtained from 
the experimental works of Yoshitake et al. (Yoshitake et al., 2012) and Zhao et al. (Zhao 


et al., 2014). The development of Tmax for cement paste is acquired from the experimental 
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work of Barde et al. (Barde et al., 2005). While it is known that the progress in hydration 
of cement paste may lead to alteration of Poisson’s ratio o”, Gy. and Gf, parameters with 
respect to time, we explicitly consider these parameters remain unchanged our numerical 
analysis indicates that the evolution of v” has a negligible effect on vm and Em development. 
Since the relative size of the fracture process zone, lez is relatively larger than the size of 
mortar RVEs (e.g., mortar RVE with characteristic length /m<15 mm), the mode-I failure 
is dominated by Thax rather than GË. (Z. P. Bazant & Planas, 1997; S. S. P. S. Shah et al., 
1995). Based on previous mechanical analysis of the early-age cement paste, the two 
parameters v = 0.18 (Kupfer & Hilsdorf, 1969), and Gr = 4144 (J/m2) (Hoover & Ulm, 
2015) are used in this work. Therefore, the damage model accounts for the time-dependent 
behavior of T? „y to address its impact on damage initiation and growth, while GÈ. is kept 
constant at 410 (J/m2). Previous mechanical analysis of the early-age cement paste found 
the development of GÈ. to be in the range of 410+40 (J/m2) (Hoover & Ulm, 2015). 

e The elastic and cohesive properties of aggregate used in the experiment (Indiana limestone) 
phase is obtained from our previous experimental study as follows (Hadi Shagerdi 
Esmaeeli, 2015): Æ” = 31.07 + 3.87 MPa, v” = 0.15, Tmax = 3.8 + 0.748 MPa, and i = 
43.4 + 16.9 (J/m’). 

e The shear components of cohesive law capture an accurate nonlinear pre-peak (dominated 
by crack initiation) and post-peak behavior (dominated by crack coalescence and 
corresponding shear failure development). To this end, the intrinsic mechanisms of opening 
and shear modes are assumed to yield to two relations of Gne ~ 10.Grc (López et al., 2008a), 
and Tmax ~ 1.4.Tmax (Eldin & Senouci, 1993). 

e The development of on and 6; parameters for cement paste, ITZ™ phases along with constant 
On and òs parameters for the fine aggregate phase are determined using the cohesive bilinear 
law (òn = 2Gtc/Tmax and òt = 2Grc/Tmax). 

It is necessary to mention that these two phases are considered to exhibit a homogenous 
brittle behavior, as the size of heterogeneities is the order of magnitudes smaller than the size of 
cement paste and fine aggregate elements; thus the cement paste and fine aggregate phases are 


treated as homogeneous and brittle materials. 
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4.3.1.1 Effect of T? „, development at early ages 


For evaluation of the existence of a unique RVE for mortar (softening material) via 
continuous and discontinuous homogenization schemes, we study the mechanical response of five 
sizes of mortar RVEs to evaluate the resulting mechanical properties along with crack patterns 
under uniaxial tension at two ages: early age of t = 16 h and a later age of t = 96 h. The RVEs has 
a 45% volume fraction of fine aggregates of a maximum size of 4.75 mm. Five RVEs of size In = 
6mm, 8 mm, 10 mm, 12 mm, and 15mm are investigated. The mortar RVEs are tested under static 
load, under a displacement control and a plane strain condition. The domain of mortar RVE is 
meshed with triangular continuum elements with quadrilateral interface elements, with element 
size as small as 0.01 mm inserted between fine aggregates. The convergence of averaged o~e 
response by increasing the size of the RVE verifies the existence of an RVE. Average strain (£) 1s 
calculated through dividing the prescribed displacement by the size of RVE, while the average 
stress (o) is similarly obtained through dividing the total reaction force (at the bottom edge of RVE) 
by the size of RVE. 

Figure 4-6 plots the o-€ curves and corresponding homogenized Ta-un curves, as well as 
localized strain patterns (crack pattern) at the final stages of loading obtained for five sizes at an 
early age. Using the continuous homogenization, the homogenized response becomes more brittle 
for larger specimens, as the averaged stresses and strains are calculated over the entire domain of 
the RVE, shown in Figure 4-6(a). As a result, the RVE size is scaled with the region outside the 
localization band. However, the slope of the linear region of pre-peak response (£”), and the peak 
response itself (Thay) make convergence as the RVE size increases. In order to evaluate the effect 
of vP development on homogenized E” and v” at t = 16 h, two values of v? = 0.15 (based on 
experimental measurement at very early ages (Byfors, 1980)) and an assumed constant value of v” 
= 0.2 are used. The results indicate that the differences remain less than 3 %, implying that a 
constant value of v” = 0.2 can be considered at a very early age. It is necessary to note that the 
solution of lm = 6 mm specimen diverges considerably from the solutions of the larger specimens, 
and thus we conclude that this specimen is too small to represent the RVE statistically. Using a 
discontinuous homogenization scheme, Figure 4-6(b) shows that the obtained homogenized 
cohesive laws are independent of the size of the RVE. The averaging is performed only over the 
localization band (Ly), which essentially takes a constant width. Figure 4-6(c) depicts the crack 


pattern in the five RVEs. We observe that the cracks have more tendency to propagate around fine 
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aggregates and across the cement paste and ITZ™ leading to complex crack patterns with more 


crack zig-zag mechanism. This is due to the fact of Tp ay < i. at t= 16 h which corresponds to 


an age of cement paste before the knee point. 
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Figure 4-6. Existence of a unique RVE for heterogeneous mortar RVEs at an early age of t = 16 h: 
(a) altering homogenized o-€ relations obtained with continuous homogenization scheme versus 
(b) homogenized Tn-un relations obtained with discontinuous homogenization scheme, and (c) 
localized strain patterns. 


To investigate the effect of mesh size on the mechanical response of mortar at t = 16 h, the 
range of mesh size studied goes from 0.1 mm to 0.3 mm. In the analyzed range, it is observed that 
the error in E” is less than 5%, in Ty", is less than 6%, and in Gj? is less than 5%. For this range 
of mesh sizes, we observe no changes in the crack pattern and number of fractured aggregates. 


Finally, we adopted an optimum mesh size of 0.15 mm that enables the best compromise between 
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the mechanical properties and the computational time. The current model has 10707 nodes with 
two degrees of freedom per node. For an explicit scheme, the time step is chosen to be 10° secs to 
satisfy the stability criterion. The computations are carried out to a total of 0.316 mm in 
displacement (until loads drop due to failure) with a total number of time steps of 326,000. 

Figure 4-7 analyzes the mechanical response of mortar RVEs at a later age of 96 h, 
including the homogenized o~e and Tn-un responses as well as the crack pattern at the final stages 
of loading. The averaged o-é curves for five RVE sizes are given in Figure 4-7(a), while the 
homogenized T;,-Un curves are shown in Figure 4-7(b). Similar to the analysis of mortar at an early 
age, the result depicts that the linear elastic responses are independent of the RVE size using 
continuous homogenization scheme (since the elastic slope and maximum tensile stress in the o-€ 
curves is E” and Tmax in one dimension, respectively). The comparison of difference of results 
between the experimental measurement of v” = 0.25 (Byfors, 1980) and an assumed constant value 
of v” = 0.2 indicates that the differences in terms of E” and v” remain less than 5 %, implying that 
a constant value of v’ = 0.2 can be considered at later age of t = 96 h. However, the softening 
responses are not objective w.r.t the RVE size. Deriving Tn-un curves independent of the RVE size 
consecutively implies that a unique RVE exists for mortar using the discontinuous homogenization 
scheme, shown in Figure 4-7(b). Upon the convergence of Tn-un curves, we quantify a unique value 
for Gir (i.e., the area under this curve) that leads to characterization of a homogenized cohesive 
law for mortar. Figure 4-7(c) depicts the crack pattern in the RVEs of five sizes. We observe that 
the cracking pattern takes straighter paths through aggregate particles with the less zig-zag 
mechanism as Tha, > (eee at t = 96 h which corresponds to an age of cement paste before the 
knee point. 

The change in cracking patterns in early-age mortar highlights the competition between 
two mechanisms, which is an important finding in this study. Comparing the values of T ay; 
TITZ" (=T? ax), and T/2,, the results indicated that the increase of T,?,,, and T/ZZ” leads to a 
significant increase in load capacity of mortar as the hydration of cement paste progresses; 
however, the ratio of two tensile strengths (TË ax/ D primarily governs the crack propagation. 
Knowing the fact that the ITZ™ phase is weak relative to aggregates and cement paste, this phase 
becomes preferred crack paths. At very early age t = 16 h, the crack impinging on the ITZ™ is likely 


to be deflected into the cement paste matrix as the condition for crack propagation in the cement 
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paste will be met at a lower load than for penetration across the fine aggregates (i.e., TE qx/ TIE, < 
1). Conversely, at a later age t = 96, the crack will tend to penetrate through fine aggregates when 


the inequality is reversed. In this study, we assess the role of Tmax considering constant values of 


GP, = GITZ™ = 410 (J/m?) and GJ? = 43.4 (J /m?). 


L, = 6mm 

l, = 8 mm 

l, =10 mm 

n 

l, = 12mm 
PONERSE "=n l = ]S mm 


l, = 6mm 
l, = 8mm 
L = ]0}0 mm 
l, = 12mm 
E a lS ie 


6, (MPa) 
o, (MPa) 


T = 





0 i 

0 0.0005 0.001 0.0015 0.002 0 0.002 0.004 0.006 0.008 0.01 0.012 
€, (mm/mm) u (mm) 
(a) (b) 





_sT A 





lmn = 6mm lm = 8mm 
l, = 4.33 mm l, = 4.33 mm 


Figure 4-7. Existence of a unique RVE for heterogeneous mortar RVEs at an early age of t = 96 h: 
(a) altering homogenized o-€ relations obtained with continuous homogenization scheme versus 


(b) homogenized Tn-un relations obtained with discontinuous homogenization scheme, and (c) 
localized strain patterns. 


4.3.1.2 Determination of E™ (t), T® x(t), and Gye 


While experiments are useful to quantify the tensile stiffness (£) and strength (Tmax), they 


provide no knowledge on the onset of failure. To this end, we employ our developed numerical 
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methodology considering the mortar as a heterogeneous composite with interfaces to interpret the 


experimental results. This model allows us to look at how the knee point in the bilinear mortar 


pra 


max and the distribution of stress 


curves (point of change in slope) is related to variable T;),,,/ 
respect to the degree of hydration of cement paste, which cannot be explicitly quantified from 
experiments. Therefore, we use the experimental Æ and TË „y values for early-age cement paste 


by Zhao et al. (Zhao et al., 2014) and Barde et al. (Barde et al., 2005), respectively, along with £" 


and Tio values for aggregate (Indiana limestone used in this study) by Esmaeeli (Hadi Shagerdi 
Esmaeeli, 2015), as shown in Figure 4-8, to evaluate the role of E and Tmax of mortar constituents 
on the early-age mechanical response of mortar. These results obtained numerically for mortar are 
compared to experimental observations of the flexural behavior of cement paste and mortar beams 
(performed by Barde et al. (Barde et al., 2005)). As previously mentioned, we consider that the 
tensile strength development of ITZ™ used in this study is identical to that of the paste due to a 
lack of experimental data on the early-age tensile strength development of ITZ™. However, to 
assess any effect of the strength of ITZ™ on the overall properties of the material we carried out a 
conservative analysis considering that the tensile strength development of ITZ™ is reduced to 50% 
of the paste (as experimentally determined by Hsu et al. [67] for later ages). The numerical results 
indicate that the tensile strength of mortar is reduced by 12.5% and 7% at t= 16 h and t = 96 h, 
respectively. However, we found that the number of fracture aggregates and the general crack 
pattern did not change significantly. While this model can easily incorporate more experimental 
data as it becomes available, it does not change the main observations discussed in this paper. 
Figure 4-8(a) illustrates the variation in Æ” and E” as a function of cement paste age (1.e., 
proportional to the degree of hydration). The slopes between E” as well as E” and t resemble each 
other where the curve representing the variation of E” is shifted to higher value respect to E? curve, 
due to the presence of fine aggregates in the mortar. Note that E depends on the behavior of the 
specimen under very low levels of loading (1.e., the slope of c-e curve in the linear elastic region). 
At these low load levels, mortar does not develop significant damage, and therefore the aggregates 
have not fractured, and as such, it is assumed not to influence the E” development. Comparing the 
numerical results and experimental data, it is observed that the percentage difference becomes less 
than 10% calculated at 16 discrete ages, except at two very early ages of 10 h (25%) and 13 h 


(27%), which implies a good agreement. 
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The next step is to predict the development of Tmax as a function of cement paste age, and 


to contrast the numerical results with experimental observation performed by Barde et al. (Barde 


et al., 2005). Given Tf ay and i... a bilinear relationship is observed for development of Thaxas 
a function of cement paste age while paste shows a linear response, shown in Figure 4-8(b). We 
observe that the general trend between the strength and age is similar at early ages (1.e., before 25 
hours) for paste and mortar (1.e., (a / EEA < 1). However, this trend significantly deviates for 
mortar once T? , (and T/7Z”) exceeds the T/S, G.e., TPax/TfS, > 1), highlighting the 
considerable role of aggregate on prediction of strength development of mortar under mode-I 
loading. Figure 4-8(b) shows that the comparison shows that the deviation between numerical and 
experimental results is less than 5%, except at two very early ages of 10 hours and 16 hours, where 
the deviation is 47% and 8.25%. With discontinuous homogenization scheme, the fracture energy 
associated with the amount of energy dissipation over the localized band has been estimated to be 
approximately Gj? =15.7 + 3.4 (J/m’). Comparing the elastic and cohesive properties quantified 
with the numerical model and experimental data, it is concluded that the employed continuous and 


discontinuous homogenization schemes can reasonably predict the mechanical properties of 


mesoscopic mortar. 
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Figure 4-8. (a) Numerical and experimental results for early-age mechanical properties 
development of mortar in terms of (a) E”, and (b) ff". 
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4.3.2 Step 2: Toward concrete composite at Level III 


This study considers the concrete material as a three-phase composite that contains coarse 
aggregate particles, surrounded by a transition zone (ITZ*) phase, and embedded in the 
homogenized Level II mortar. Similar to mortar, we consider a zero width for ITZ in concrete RVE 
to avoid the expensive computational cost, as the average width of ITZ is 1/40" of the average 
element size. Figure 4-8(a) indicates that the tensile stiffness of early-age mortar is most likely to 
approach the tensile stiffness of coarse aggregate as the hydration of cement paste progresses at an 
early age (i.e., E™ ~ E°*). Constantinides et al. (Constantinides & Ulm, 2004) suggested that the 
key morphological parameters of inclusions including the size gradation, spatial distribution, and 
specific fraction, play important roles in the determination of local stress distribution, micro-crack 
initiation and coalescence resulting in the macroscopic failure of concrete (Wriggers & Moftah, 
2006). Hence, the initial step to developing a suitable input for this mesomechanical model is to 
acquire appropriate parameters mentioned above for coarse aggregates within the concrete RVE. 
To this end, 3D mesostructures of concrete with randomly distributed coarse aggregates of 
spherical shape are generated, where the fraction and size distribution of aggregates is controlled. 

Prior experimental (F. H. Wittmann et al., 1988) and numerical studies (Lopez et al., 2008a) 
have shown that volume fraction has the most pronounced influence on the mode-I fracture 
behavior such as load capacity (1.e., tensile strength), cracking tendency, or fracture energy. Here, 
we develop an in-house composite material generator program composed of spherical aggregate 
generation and packing where the specific size distribution and the specific fraction of aggregates 
(fea = 38%) are determined experimentally by Barde et al. (Barde et al., 2005) based on the mixture 
proportioning for concrete. To place a spherical aggregate at an arbitrary position within the 
mesostructure, three requirements should be met: (1) the whole sphere needs to be placed within 
the boundaries of the mesostructure with no intersection; (2) no overlap with previously placed 
spheres should occur; and (3) A limitation on two minimum distances between the edge of a sphere 
and the boundaries of the mesostructure, as well as the gap width between two adjacent spheres 
should be set. Figure 4-9 depicts the flowchart of processes of packing of aggregate, including 
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“input,” “taking,” and “placing” (X. Wang et al., 2015). The “input” process examines the 
prominent parameters, including specimen size, aggregate size gradation, and aggregate volume 
fraction, for generating an RVE with randomly distributed aggregates. According to the random 


size description, the “taking” process produces a single aggregate particle. Next, the aggregate 
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placement into the predefined location in an arbitrary fashion 1s determined through the “placing” 
process, in which prescribed physical constraints (intersection and overlap checking) are 
implemented. Three following conditions need to be considered to check even if there is any 
intersect or overlap among randomly distributed aggregates: (1) the entire domain of aggregate 
particle must be fully embedded within the box, (2) no overlap should exist with previously 
embedded aggregate particles, and (3) the minimum gap width between the surface of an aggregate 
and the edge of the box as well as the minimum distance between two neighboring aggregate 


particles should be set. 


Read geometrical parameters: 
Specimen size, aggregate volume 


fraction, aggregate gradation 


i“ aggregate size range: 
Generate an aggregate with random 
size and location 


| 
Yes 
Intersection? 
No 
satisfied? 
Yes 


Figure 4-9. Flowchart of RVE generation with random aggregate distribution, adapted from 
(X. Wang et al., 2015). 
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4.3.2.1 The representative fraction of aggregates 


Afterward, 2D planar slices are extracted from 3D mesostructures to obtain a specific 
fraction of aggregates in 2D mesostructures statistically. The statistical analysis indicates that a 
skewed normal distribution can well describe the fraction distribution of aggregates with a mean 
value fea = 45% and a standard deviation of 3% for a 2D mesostructure. Afterward, 2D 
mesostructures with the specific fraction of coarse aggregates acquired from statistical analysis are 
extracted from the 3D mesostructures of concrete with randomly distributed coarse aggregates of 
polyhedral (realistic) shape. The newly generated 2D and 3D mesostructures are tested under 
uniaxial tension to assess the validity of numerical results by contrasting the mechanical response 
of the 2D model with the 3D model, along with experimental and analytical results. The need for 
Statistical analysis is since assumption of an identical fraction of coarse aggregates for 2D and 3D 
specimens is a crude approximation as the 2D model does not capture the constraint effect in the 
thickness direction that may bring significant underestimation of pre-peak and post-peak responses. 
We also assume that the failure of mortar and coarse aggregate phases obey Weibull statistics. It 
means that they contain uniformly distributed defects using Weibull distribution in which their 
strength 1s determined based on the number of defects subjected to tensile stresses with high values. 

Prior experimental and numerical studies have shown that volume fraction has the most 
pronounced influence on the mode-I fracture behavior such as load capacity (1.e., tensile strength), 
cracking tendency, or fracture energy. Here, we develop an in-house composite material generator 
program composed of spherical aggregate generation and packing where the specific size 
distribution and the specific fraction of aggregates (fea = 38%) are determined experimentally 
based on the mixture proportioning for concrete. To place a spherical aggregate at an arbitrary 
position within the mesostructure, three requirements should be met: (1) the whole sphere needs 
to be placed within the boundaries of the mesostructure with no intersection; (2) no overlap with 
previously placed spheres should occur; and (3) A limitation on two minimum distances between 
the edge of a sphere and the boundaries of the mesostructure, as well as the gap width between 
two adjacent spheres should be set. More information on a comprehensive procedure for the 
generation of 3D mesostructure of concrete with spatially distributed spherical aggregates. Based 
on the developed procedure, a series of 3D mesostructures (50 x 50 x 50 mm”) consisting of mortar 
and spherical aggregates are generated and shown in Figure 4-10(a). The inset in this figure 


displays a typical size distribution of coarse aggregates obtained from INDOT specification [67]. 
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Afterward, 2D planar slices from each one of 100 constructed mesostructures are extracted every 
1 mm along with x, y, and z axes. Note that the first slices extracted from each side of the 3D 
mesostructure are neglected in our analysis due to the limited space between the edge of the coarse 
aggregate and the boundary of the mesostructure (1.e., a total number of 14100 2D slices are 
extracted). The fraction distribution of aggregates can be well described by a skewed normal 
distribution with a mean value of 45% and a standard deviation of 3% for a 2D mesostructure), as 


shown in Figure 4-10(b). 
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Figure 4-10. (a) 3D mesostructure of concrete consisting of mortar and spherical aggregates with 
a specific fraction fea = 38% spatially distributed, (b) statistical analysis of the highest probability 
of occurrence of aggregate fraction in a 2D mesostructure extracted from 3D mesostructure of 
concrete. It is observed that a skewed normal distribution (fitted red curve) describes the 
distribution of aggregate fraction with a mean value fea = 45% in 2D mesostructures extracted from 
3D mesostructure. 
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4.3.2.2 Generation of 2D and 3D concrete mesostructures 


To investigate the validity of the proposed representative area fraction of aggregate, the 
mechanical response of 2D and 3D concrete RVEs under pure tension are assessed. To this end, a 
series of cubic 3D concrete RVEs (50 x 50 x 50 mm”) containing arbitrarily distributed coarse 
aggregates of realistic (polyhedron) shape with a volume fraction of 38% and maximum aggregate 
size of 25 mm is generated, shown in Figure 4-11. The generated convex polyhedrons are made 
via the digitization of real crushed coarse aggregates with nominal sizes between 5 mm and 25 mm 
(in accordance with ASTM C-33 [68]). The main idea is to generate polyhedron aggregates with 
arbitrary orientations in a random arrangement until the target fraction of aggregates is reached. 
For this case, the geometrically arranged aggregates being placed within the concrete RVE need 
to satisfy two conditions: (1) despite the procedure of placement of spherical aggregates, there is 
no distance between the surface of a polyhedral aggregate and the edge of the concrete specimen 
to avoid the decrease of aggregate volume fraction on the boundaries; (2) if an aggregate crosses 
through the boundary of concrete RVE, the volume of the portion of aggregate being placed will 
be considered for calculation of aggregate fraction; (3) the minimum distance between any two 
aggregates are considered to be 0.2 mm. Figure 4-11(a) shows an example of 3D concrete RVE. 
Similar to the previous section, 2D planar slices with aggregate area fractions of 45%, shown in 
Figure 4-11(b) and (c), and 38%, as shown in Figure 4-11(d) and (e), are extracted from the 3D 


RVE to assess the validity of proposed representative aggregate area fraction. 
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Figure 4-11. (a) A three-dimensional RVE of concrete containing irregular shape of coarse 
aggregates with fca = 38%; two-dimensional slices extracted from three-dimensional concrete RVE 
containing aggregates with (b, c) fea = 38%; (d, e) fea = 45%. 


To account for the effect of pre-existing cracks in the mortar, coarse aggregate, and ITZ* 
phases (these materials contain voids with dimensions up to a few millimeters) on the resulting 
peak stress, a Weibull distribution of strength for concrete constituents is employed [69-71]. 
Theoretically, this statistical approach enables the definition of strength distribution of concrete 
material considering a large variability in tensile strength of concrete constituents (mortar and 
aggregate) due to the variability in the distribution of crack size, shape, and orientation concerning 
the loading direction. It should be noted that ductile materials exhibit a very narrow strength 
distribution curve, which can be represented by Gaussian or normal distribution. In contrast, brittle 
and quasi-brittle materials exhibit a fairly broad strength distribution which can be described by 
Weibull distribution. This study considers a spatial distribution of strength, Tmax for interface 
elements as fracture energy per area, Gre is assumed to remain unchanged. Zhou et al. (F. Zhou & 
Molinari, 2004) described that the strength of the interface element (Tinax) strongly depends on the 


local length scale, namely an interface element volume (V). Therefore, Weibull distribution 
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enables the characterization of the strength of mortar and ITZ* from the numerical results, as well 
as the strength of coarse aggregate from the flexural experiment, and the use of this information to 
predict the peak stress and the onset of failure in the mesostructure of concrete. Considering a 
volume V consisted of equivalent reference volume Vo and a reference tensile strength fs the 
strength distribution for interface elements within concrete constituents is described by equation 
4-6: 


f Tmax) = an exp (- = = ) 4-6 





where m is the Weibull modulus that controls the shape of the distribution, the exponent m specifies 
how the strength is distributed around the average strength Tmax. The shape parameter generally 
takes on a value between 5 and 20 for the case of brittle materials, including mortar and aggregate 
materials. Within our model, m = 20 describes a practical defect-free material, while m = 5 
determines a thoroughly heterogeneous material. We have two values of m = 5 and m = 7 for 
mortar (ITZ*) and aggregate phases to capture the influence of two various defect distributions on 
the resulting strength and crack patterns. 

Equation 4-6 indicates that an interface element with a smaller V compared to Vo exhibits 
a lower probability of failing subjected to given tensile stress as it is less probable to have defects 


1/m 
l V . i 
(een Tna X Je (—) ). In order to obtain the material parameters used in equation 4-6, we first 
0 


assess the influence of length scale on Tmax of mortar, ITZ‘, and coarse aggregate in the concrete 
mesostructure. Because Tmax of mortar is acquired from numerical simulation of Level IJ mortar 
RVE with an average interface element size of 1.5 mm, then we adopt Vo = 1.5 mm? (simulations 


are performed under plane strain condition). Considering an average finite element size in 2D 


Cc 
mesostructure of concrete to be 5 mm (i.e., V = 5 mm’) and m = 5, es is set to be equal to 


1/m 

EZ V A I 

i as (—) ~ 1. For coarse aggregate, the higher values of tensile stress occur over a 
0 


region in which the typical size of material inhomogeneities, known as, fracture process zone (FPZ) 
is characterized. Bazant et al. (Z. P. Bazant & Planas, 1997) suggested that large specimens tend 


to fail in a brittle manner when FPZ size << specimen dimension. Therefore, we simply assume 


1/m 
that the problem of size effect for coarse aggregate phase can be neglected (1.e., (—) ~ 1)as 
0 
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the length of FPZ (= 38 mm) 1s significantly smaller than the specimen size tested in a four-point 
flexural experiment (152.4 x 152.4 x 457.2 mm’). 

The effect of m parameter, for instance, on the spatial distribution of strength within concrete 
mesostructure, is shown in Figure 4-12. The left and right insets show the strength variation and 
spatial distribution within a random region of concrete RVE for selected shape parameters m = 5 
and m = 7 for mortar and aggregate phases, respectively (blue color associates with lower strength 


values, and red color associates with higher strength values). 
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Figure 4-12. A plot of two Weibull distributions of tensile strength with shape parameters 
of m = 5 and m = 7 within concrete RVE. Histogram for Tmax, the color map presents the 
ratio of the tensile strength of interface elements respect to constant tensile strength, Tmax. 


The early-age elastic properties of the three phases in concrete and the cohesive properties of 
the interface elements are defined as follows: 

e Regarding the Level III concrete (and ITZ‘), the development of E™(t) and Thax (t) as well 

as constant Poisson’s ratio, v™ and Gj? are obtained from the homogenization of the 


mechanical response of Level II mortar. The elastic and cohesive properties of coarse 
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aggregate are identical to the properties determined for fine aggregate, as shown in Table 
1. Similarly, the relations for shear components of cohesive law for concrete constituents 
are identical to the relations defined for mortar constituents. 

e The development of on and ò parameters for mortar, ITZ* phases along with constant òn and 
or parameters for the coarse aggregate phase are determined using the bilinear cohesive law 


(On = 2Grīc/ Tros and Òt = 2GiItc/ Tmax). 


4.3.2.3 Effect of the aggregate fraction on the mechanical response of concrete 
mesostructure 


To evaluate the effect of fea on the tensile mechanical response of concrete including E“ and 
Tmax at a very early age, t = 16 h, numerical results for 2D mesostructures with statistically 
representative fea = 45%, and 2D/3D mesostructures with fea = 38% are analyzed under uniaxial 
tension and contrasted with analytical and experimental solutions, shown in Figure 4-13. Three 
mesostructures of 2D concrete with fea = 38%, 2D concrete with fea = 45%, and 3D concrete with 


fea = 38% are investigated. The material parameters of the mesostructural constituents can be found 


in Table 4-2. 


Table 4-2. Material parameters of the continuum and interface elements within concrete 
RVE. 


Mechanical t= 16h t= 96h 
Unt aaaaaaaaaaaaaaamammmmmmmmmmaħįÃ ca 
Property m ITZ m ITZ 
E [GPa] 14.02 - 36.81 - 31.074 
V [-] 0.18 - 0.18 - 0.18 
Gr [J/m"] 15.7 15.7 7 15.7 43.4 
Gi [J/m?] 157.0 157.0 157.0 157.0 434 
Tmax [MPa] 3.24 3.24 5.862 5.862 3.8 
Tmax [MPa] 4.53 4.53 8.8 8.8 5T 
Ôn [mm] 0.0096 0.0096 0.0053 0.0053 0.0013 
Ôt [mm] 0.069 0.069 0.035 0.035 0.015 
lez [mm] 20.96 20.96 16.82 16.82 38.103 
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While numerical simulations widely quantify tensile stiffness of heterogeneous composites 
containing inclusions with complex geometries, two analytical solutions of Mori-Tanaka (MT) 
and self-consistent method (SCM) were proposed based on mean-field homogenization schemes 
(R. Christensen et al., 1992; Hori & Nemat-Nasser, 1994). These two analytical methods were 
established based on the linear elasticity theory that estimates the effective tensile stiffness of 
composite using elasticity solution of Eshelby’s for inclusions placed in an infinite matrix. These 
two models are best suited for the given homogenization problem, as their predictive capabilities 
are effective when (1) the composite contains a high fraction of inclusions, and (2) the matrix has 
a relatively low stiffness compared to inclusions (R. M. Christensen, 1990). This allows us to 
compare E“ of concrete mesostructures with that of idealized elastic two-phase solid, showing how 
the numerical model not only gives logical results but is also capable of fitting analytical solution 
and proper experimental data (Antico et al., 2015), as shown in Figure 4-13(a). It is observed that 
the mean value of E“ predicted by 2D and 3D mesostructures with fea = 38% are 16.2 and 19.1 
GPa, respectively (1.e., the relatively more realistic 3D specimen predicts 15.1 % higher E* value). 
This significant difference can be alleviated to a relatively negligible difference of less than 5% 
by increasing fca of 2D mesostructures to 45%. Also, the simulated E“ of 2D mesostructures with 
fea = 45% 1s following MT and SCM solutions for a heterogeneous composite with an inclusion 
volume fraction of 45%, as well as experimental measurements. 

Next, we evaluate the capability of the 2D mesostructure with fea = 45% on the prediction 
of Thax aS well as damage mechanisms. We then carried out numerical simulations with two 
values of m varied between 5 and 7. The mean peak stresses from 2D and 3D modeling against fea 
are presented in Figure 4-13(b). The mean peak stresses obtained from the analysis of 2D and 3D 
mesostructures with fea = 38% are 1.57 and 1.65 MPa, respectively, 1.e., the relatively more 
realistic 3D mesostructure estimates 4.8 % higher peak stress. It is concluded that the rise in the 
peak stress in the 3D model compared to the 2D model with identical fca is associated with the 
constraint effect through-thickness direction considered in 3D modeling, not accounted by the 2D 
model. Liu et al. (L. Liu et al., 2016) also confirmed that this effect yields a uniform distribution 
of pre-peak meso-cracks in 2D, as the spatial distribution of aggregates may have a lesser influence 
on the pre-peak response in the 3D model. To alleviate the 2D effect, we analyze the peak response 
of 2D mesostructures with fea = 45% (rather than fea = 38%), where the difference between the 


mean peak stresses of 2D and 3D models reduces to less than 4 %. One typical macro-crack pattern 
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with only one dominant crack propagated across ITZ* and mortar phases is observed from 2D and 
3D simulations and illustrated in the inset of Figure 4-13(b). Afterward, the peak stresses obtained 
from 2D mesostructure with fea = 45% and 3D mesostructure with fea = 38% are contrasted with 
experimental observation (Barde et al., 2005). The experimental work investigated the mechanical 
response of concrete beams subjected to a four-point flexural loading, where a constant bending 
moment is induced between the inner two loadings with no effect of shear stress. As a result, the 
local stress and strain fields on the bottom surface between two loadings are subjected to equivalent 
tensile stress which implies that failure mechanism is dominated by mode-I once maximum 
principal stress exceeds Thay. Compared to experimental observations, the 2D and 3D modellings 
demonstrate a prediction of mean peak stress with a difference of 2% and 6%, respectively, 1.e. a 


better agreement is suggested between 2D model with fea = 45% and experiment. 
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Figure 4-13. (a) Analytical, numerical and experimental results for homogenized E* parameter of 
concrete at t = 16h. A 2D model with fea = 45% produces a negligible uncertainty of 5% in 
prediction of E* parameter contrasted with the experimental measurement; (b) Comparison of 
numerical and experimental results for homogenized f; parameter of concrete at t = 16 h. Similar 
to the prediction of the E parameter, the difference between the mean peak stress predicted by 2D 
with fea = 45% (rather than 38%) becomes less than 4 %. 
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4.3.2.4 Development of mechanical properties of concrete 


In this section, the mechanical contribution of constituents on mechanical properties 
development of composite concrete 1s analyzed. To this end, we extend the numerical results, 
illustrated in Figure 4-13, to include the aging effect of mortar along with the mechanical 
contribution of coarse aggregate using the same continuous homogenization scheme. To verify the 
validity of the numerical model, the development of mechanical properties of concrete including 
E° and ff are quantified at early ages (16 discrete ages), and then contrasted with proper 
experimental observations (Antico et al., 2015) and (Barde et al., 2005), respectively. Following 
the experimental protocol (Barde et al., 2005), the four-point flexural loading is simulated using 4 
rigid rollers having frictionless contact between the rollers and the concrete beam (508 x 152.4 
mm’), displayed in Figure 4-14(a). The images from the cross-section of the concrete beam tested 
in the experiment are used to obtain the real 2D RVEs of concrete using a computed tomography 
approach (Hadi Shagerdi Esmaeeli, 2015). To reduce the computational expense, only one RVE is 
embedded on the midspan of the beam, which consists of continuum elements of different phases 
(e.g., mortar and aggregate) along with interfacial elements inserted between them. For the 
elements outside the RVE, continuum elements with homogenized elastic properties of E and v 
are obtained from the response of RVEs with no interface elements tested under uniaxial tension. 
The elastic properties of the continuum elements are obtained from performing separate 
simulations of concrete RVE, where the development of E“ is illustrated in Figure 4-14(a) and the 
constant Poisson’s ratio v° is determined to be equal to 0.17. It is also observed that the micro- 
crack coalescence begins from the bottom edge approaching to the top edge of RVE, while the 
resulting response is still in the pre-peak region. In order to avoid this boundary effect on the 
fracture process, RVEs of sizes le = 60 mm, 80 mm, 100 mm, and 120 mm with fea = 45% and 
maximum aggregate size of 25 mm are investigated. The elastic and cohesive properties of mortar 
(and ITZ*) phases are acquired from homogenization of Level II mortar as well as the coarse 
aggregate properties are obtained from our previous experimental investigation (Hadi Shagerdi 
Esmaeeli, 2015). The concrete beam is loaded under displacement-controlled plane strain 
condition. 

To verify the existence of an RVE, the mechanical response of concrete beams in terms of 
load-midspan displacement (P-ò) relations (measured through the contact between rollers and 


beam located at the top edge) with increasing the size of the RVEs are investigated at an early age 
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of t= 16 h and a later age of t = 96 h. Figure 4-14(b) and Figure 4-14(c) show that macroscopic 
load-displacement curves converge as the size of the RVE increases. It should be considered that 
the 60 x 60 mm? specimen is too small to statistically represent the mesostructure as its solution 
diverges considerably from the solutions of the larger specimens. The comparison of results shows 
that the pre-peak response of the 80 x 80 mm”, 100 x 100 mm”, and 120 x 120 mm” specimens are 
not on top of each other or very closed to each other due to considering only one realization for a 


certain size. Therefore, the macroscopic pre-peak responses for these three specimens begin to 


2 


converge upon increasing the mesostructure size. Thus the 80 x 80 mm specimens will be 


considered as the representative 2D RVE for concrete in the following simulations. 
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Figure 4-14. Force-displacement curves of (a) concrete beams containing four sizes of RVE under 
four-point flexural loading (b) at an early age of t = 16 h, and (c) at a later age of t = 96 h. Sequential 
macro-crack pattern evolution of concrete RVE at (d, e, f) t= 16 h; and (g, h, 1) t= 96 A. 
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Figure 4-14 continued 
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A mesh size range from 0.5 mm to 2 mm is analyzed to investigate the effect of mesh size 


on the mechanical response of concrete at t = 16 h. The analyzed range indicates that the error in 
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E° is less than 7% and in Thax is less than 4%. For this range of mesh sizes, we observe no changes 
in the crack pattern and number of fractured aggregates. To achieve a balance between the accuracy 
of results and computational cost, we selected a mesh size of | mm for all the simulations of Level 
I concrete. The current model has 36427 nodes with two degrees of freedom per node. For an 
explicit scheme, the time step is chosen to be 10°° secs. The computations are carried out to a total 
of 0.0629 mm in displacement (until loads drop due to failure) with a total number of time steps of 
656,000. 

The detailed pictures of two typical macro-crack patterns with only one dominant crack are 
illustrated at various loading stages in Figure 4-14: (1) At early ages (.e., before 25 h), the 
numerical modeling suggests that the damage mechanism is dominated by mortar fracture, where 
the number of damaged aggregate particles is three highlighted in pink color, as sketched in Figure 
4-14(d-f). In the initial loading stages, a significant number of micro-cracks are observed to form 
at ITZ‘ (Figure 4-14(d)). As loading progresses, several micro-cracks coalescence into one 
predominant macro-crack, whereas redistribution causes the other micro-cracks to arrest (Figure 
4-14(e)). The single dominant crack progressively propagates through ITZ* and mortar elements 
indicating a rapid softening behavior (Figure 4-14(f)); (2) At later ages (1.e., after 25 h), the 
numerical modeling illustrates an alteration in fracture behavior of concrete as the damage 
mechanism is dominated by aggregate fracture, where a single dominant crack propagates across 
seven aggregate particles, as sketched in Figure 4-14(g-1). Looking at the details of mesostructural 
cracking, one can see how micro-cracks start developing in a distributed manner all over the 
specimen, especially at the aggregate interfaces. At a certain stage of the loading, some micro- 
cracks initiate at ITZ* and aggregate phases which continue opening and getting interconnected as 
the loading increases, as sketched in Figure 4-14(g) and (h). At the softening region, the 
localization process 1s spontaneously triggered, and the process leads finally to only one dominant 
crack cutting through the RVE into two big pieces. The results also exhibit other realistic features 
such as crack branching at an early age. It clearly shows that as Thx > Tmax(t) at early ages, 
cracks have more tendency to propagate across the mortar phase and do not deflect through 
aggregates that leads to a branching mechanism. This mechanism occurs even within small mortar 
regions between adjacent aggregate particles or corners. 

Here, we quantify the early-age tensile stiffness and strength at discrete ages using 


continuous homogenization scheme, and compare them with experimental observations, shown in 
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Figure 4-15(a), and Figure 4-15(b), respectively. The comparisons presented in Figure 4-15 also 
include the range of experimental values of E™ and T,;4, (highlighted in pink) obtained from (Hadi 
Shagerdi Esmaeeli, 2015). Figure 4-15(a) shows the relationship between E* and the age of 
concrete. The numerical results show that the development of E“ is similar to the development of 
E”. This similar trend of tensile stiffness development for early-age mortar and concrete specimens 
as the concrete has not developed significant damage. Therefore, the aggregates have not fractured, 
and as such, do not influence the tensile stiffness development. Comparing the numerical solution 
and experimental observation, it is seen that the percentage difference becomes less than 10% 
calculated at discrete ages, except at two ages of 19 h (32%) and 37 h (26%), which implies a 
promising agreement. Next, given the tensile strength of concrete constituents, the development 
of strength in mode-I as a function of age for concrete material obtained from numerical solution 
and experimental observations is shown in Figure 4-15(b). It can be observed that the slope 
between the strength and age significantly alters once the strength of the mortar (ITZ) phase 
exceeds the strength of the aggregate phase. The comparison between numerical and experimental 


results shows that the variation is less than 6%, which indicates an excellent agreement. 
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Figure 4-15. (a) Numerical and experimental results for mechanical properties development of 
early age concrete in terms of (a) tensile stiffness E£“, and (b) tensile strength ff and number of 
fractured aggregate particles in concrete RVE. 


It is crucial to understand the competing mechanisms associated with the energy dissipation 


in heterogeneous mesostructure of cementitious composites to address their remarkable early-age 
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mechanical properties development. Our study indicates that two competing mechanisms exist 
between energy dissipation through the failure of mortar and ITZ" phases at an early age, and the 
energy dissipated at the aggregate phase at later ages. To verify the hypothesis considered in the 
experiment, this behavior is currently considered to be attributed to (1) weak interfaces of mortar 
relative to coarse aggregate interfaces which become the preferred crack paths at early ages, when 
the development of tensile strength of mortar due to aging effect leads to a similar strength gain 
for concrete; and (2) weak interfaces for aggregate relative mortar interfaces leads to alteration in 
the preferred crack path which leads to a considerable increase of the number of fractured 
aggregates and therefore limiting the strength development governed by the aggregate strength. 
The numerical framework developed in this work could be used to predict the correct tensile (or 
flexural) strength development for cementitious composites at an early age, where the presence of 


aggregates is considered to modify the maturity-strength gain relationship. 


4.4 Conclusion 


This paper investigated the early-age tensile stiffness and strength development of 
cementitious composites such as mortar and concrete through the finite element model, which were 
compared with experimental measurements. Following the work performed by Barde (Barde et al., 
2005), the flexural results indicated that a linear relationship existed between the strength and 
degree of hydration for paste specimens. However, a bilinear relationship is observed as the 
aggregate was added to the paste (1.e., mortar and concrete). The rate of development of tensile 
stiffness contrasts with that of strength, as the fracture of aggregates does not affect the stiffness 
as a function of the degree of hydration. A two-step homogenization model multiscale method 
consisted of continuous and discontinuous homogenization schemes for multiscale modeling of 
failure of quasi-brittle mortar (Level II) and concrete (Level III) having heterogeneous 
mesostructures were used. 

e Level II mortar: Mesostructure of mortar contains fine aggregates dispersed in cement 
paste, equipped with intra-phase (fine aggregate/fine aggregate and cement paste/cement 
paste) and inter-phase (fine aggregate/cement paste) cohesive elements. While the 
geometry for mortar is not explicitly considered, the mechanical properties of mortar 
constituents strongly affect the mechanical response and crack propagation in mortar. The 


continuous homogenization scheme has been used for quantification of pre-peak and peak 
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stress response of mortar. At the same time, the behavior of a macroscopic cohesive crack 
is defined upon microscopic localization using the discontinuous homogenization scheme. 
The comparison between numerical solutions and experimental observations obtained for 
tensile stiffness and strength development at discrete early ages showed that the variations 
are less than 5%. With discontinuous homogenization scheme, the fracture energy 
associated with the amount of energy dissipation over the localized band has been 
estimated to be approximately Gj? = 15.7 + 3.4 (J/m’). Comparing the elastic and cohesive 
properties quantified with the numerical model and experimental data, it is concluded that 
the employed continuous and discontinuous homogenization schemes can perfectly predict 
the mechanical properties of mesoscopic mortar. We also observed that the pre-peak 
response of mortar is strongly sensitive to the tensile stiffness of their constituents. In 
contrast, the differences between constituents’ strengths lead to govern the micro-crack 
coalescence and macro-crack pattern. We validated the hypothesis of this work by showing 
the influence of aggregate fracture on the knee point is primarily controlled by the 
competition between the energy dissipation mechanism of the cement paste at very early 
ages and the energy dissipation mechanism of the aggregate phase at later ages. 

Level III concrete: Mesostructure of concrete contains coarse aggregates embedded in 
mortar, equipped with intra-phase cohesive elements (coarse aggregate/coarse aggregate 
and mortar/mortar), as well as and inter-phase cohesive element (coarse aggregate/mortar). 
An efficient algorithm for generating three-dimensional mesostructures with the known 
volume fraction of aggregate (fca = 38% used in (Barde et al., 2005)) of spherical shapes 
were developed to perform a statistical analysis to determine a representative area fraction 
of aggregate (fca = 45%) for two-dimensional mesostructures of concrete. Afterward, a 
three-dimensional mesostructure of concrete containing realistic polyhedral aggregates 
with fea = 38% is generated, where four two-dimensional mesostructures with fea = 38% 
and 45% are extracted to compare their pre-peak response and peak stress with analytical 
and experimental observations at an early age of t= 16 h. We observed that an increase in 
area fraction of aggregate in two-dimensional specimen alleviates the 2D effect, as the 
variation between predicted numerical results and experimental data decreased to less than 
5%. Hence, the images from the cross-section of the concrete beam tested in the experiment 


are used to obtain the real 2D RVE of concrete using a computed tomography approach. 
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With a continuous homogenization scheme, the uncertainty of results obtained from the 
predictive model for early-age tensile stiffness and tensile strength development was less 
than 10% and 6%, respectively, that indicates an outstanding agreement. Similar to mortar, 
we observed that by understanding the composite nature of early age concrete, tensile 
stiffness and strength development were captured. It was observed that the tensile stiffness 
of mortar and concrete increases at a very similar rate, while the presence of coarse 
aggregates influences the rate of strength development. It was shown that the hydration of 
the paste phase describes the strength gain at an early age; however, it is primarily governed 
by aggregates after a certain degree of hydration is reached. Aggregates do not influence 
the tensile stiffness development similar to how they influence the flexural strength 


development. 
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5. TOWARDS FABRICATION OF ENHANCED DAMAGE TOLERANT 
ARCHITECTURED CEMENT-BASED MATERIAL 


This chapter contains experimental data and numerical simulations collected by Hadi S. 
Esmaeeli. Data conceptualization was a collaborative effort with Pablo D. Zavattier1. The 
manuscript was written by Hadi S. Esmaeeli and edited by Pablo D. Zavattieri. 


5.1 Problem definition 


A recent development in additive manufacturing (AM) of cementitious composites have 
demonstrated that this technology has the potential to address some key challenges such as 
fabrication of multifunctional structures with enhanced controllability over architecture, 
composition, and incorporation of components (e.g., reinforcement, sensors, and nanomaterials) 
in the construction industry compared with those resulting from conventional casting techniques. 
This technology also opens new avenues for overcoming durability issues by fabricating a 
carefully tailored heterogeneous structure at an architectural scale. A primary source of premature 
cracking of cementitious composites is restrained shrinkage with early-onset, due to temperature 
variation, self-desiccation, and non-uniform drying. In this chapter, we assess the mechanical 
performance of 3D-printed cementitious composites along with the role and characteristics of 
processing-induced heterogeneous interfaces in controlling and mitigating shrinkage cracking. 
Particularly, this work focuses on utilizing weak interfaces between filaments to promote a unique 
damage localization mechanism, where the damage propagation to neighboring layers is inhibited. 
For further performance improvement, we employ a novel architecture inspired by Bouligand 
design to promote a damage delocalization mechanism at unrestrained layers. This novel 
integration of two damage mechanisms leads to enhancing the damage tolerant and fracture 
characteristics, which allows the cementitious composites with early-age autogenous and drying 
shrinkage behavior to accommodate the damage evolution, without sacrificing the inherent early- 


age strength development. 


5.2 Literature review 


Ceramic-based materials, such as traditionally cast cementitious materials, fall into the 


category of material classes with the highest strength and stiffness (Ashby & Cebon, 1993). A 
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desirable requirement for the majority of structural materials is that they exhibit attainment of high 
strength and toughness properties simultaneously (generally, these two specific material properties 
are mutually exclusive (Robert O Ritchie, 2011)). The strength of material represents its invariable 
resistance characteristic when a material undergoes non-recoverable (inelastic or plastic) 
deformation, and toughness represents the invariable energy required to cause fracture the material. 
Owing to strong and directional bonding between constitutive phases at atom-level in ceramics, 
these materials are unable to locally dissipate high stresses developed when they undergo limited 
unrecoverable deformation. Conversely, materials with lower strength exhibit higher deformation 
and tend to be tougher (shown in Figure 5-1). 

While it has been generally thought that ceramics with iono-covalent atomic bonding is 
unlikely to undergo plasticity corresponding to the dislocation movements in crystalline materials, 
prior research has discovered numerous types of damage-tolerant biological materials, such as 
glass sponges, tooth enamel, and nacreous part of seashells, which their primary constituents 
(hydroxyapatite, aragonite, silica glass) inherent brittle behavior surprisingly (Dunlop & Fratzl, 
2010; Y. F. Gao & Bower, 2004; Ji & Gao, 2010; Launey et al., 2010). This counter-intuitive 
behavior has led to many investigations performed on biomaterials (Francois Barthelat et al., 2016; 
Chen et al., 2018; Eder et al., 2018; Faber et al., 2018; Fratzl & Weinkamer, 2007; Grunenfelder 
et al., 2018; Kokkinis et al., 2015; S. Ling et al., 2018; Zengqian Liu et al., 2017; Meyers et al., 
2013; Naleway et al., 2016; C. Ortiz & Boyce, 2008; Su et al., 2016; Weaver et al., 2012; W. Yang 
et al., 2019; Yaraghi & Kisailus, 2018; Z. Yin et al., 2019), where the dissipating mechanisms fall 
into two distinct categories: 1) the increase of crack growth resistance occurs at the traveling crack 
tip occurring within the fracture process zone (an inherent material property) due to inelastic 
deformation of soft constituents, micro-cracking, and fiber bridging and pullout (F Barthelat & 
Espinosa, 2007; Espinosa et al., 2009; Launey & Ritchie, 2009; Rabiei et al., 2010); 2) crack 
growth resistance increases by extending the crack size due to fracture process zone lengthening, 
or microcrack deflection through weak interface network (Launey et al., 2010; Launey & Ritchie, 
2009; R O Ritchie, 1988). 

For most of the anisotropic bioinspired designs, the key properties of interfacial materials 
are different than of constituents which is common among various biological systems: 1) 
interfacial regions are exhibit more compliant behavior and therefore are more likely to deform 


than the constituent, and 2) interfaces are typically weaker and attain lower flaw-resistant 
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properties (Zengqian Liu et al., 2020; Moini et al., 2018). These physical characteristics are solid 
foundations for promoting a unique toughness mechanism by deflecting and twisting the crack 
through the complexity of (micro-)structural orientations (Currey & Kohn, 1976; Kamat et al., 
2000; Mayer, 2005; Suksangpanya et al., 2018). Therefore, there 1s a dramatic rise of interest in 
utilizing additive manufacturing (AM) and multiscale hierarchal design of multifunctional and 
multilayered composites with unique flaw-tolerant and damage-resistant properties (Compton & 
Lewis, 2014; Duoss et al., 2014; Gosselin et al., 2016; Khoshnevis et al., 2001; Wegst et al., 2015). 
Given the flexibility of AM technique on combining arbitrary hierarchal structural orientation and 
interfacial materials of distinct properties, current research on adaptation of AM technology to the 
production of cement-based materials focus on eliminating the weak interfacial regions to improve 
their resulting mechanical properties and durability performance (Khan et al., 2020; Nerella et al., 
2020; Panda et al., 2017; Zareiyan & Khoshnevis, 2017). Unlike the approach of elimination of 
AM-induced interfacial regions in cement-based materials, we surmise that the fabrication of 
weaker interfaces combined with an optimized structural orientation (inspired by nature) 
encompassing anisotropy leads to desirable damage-resistant properties under intrinsic shrinkage 
loading. We hypothesize that harnessing the unique physical properties of heterogeneous 
interfaces (more compliant and weaker when compared to filaments) introduces the easy 
propagation of cracks along with interfaces and crack growth inhibition into unrestrained 
neighboring layers for promoting the flaw-resistance properties, without sacrificing the strength of 


the cementitious composite. 
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Figure 5-1. Fracture toughness-strength relation of engineering materials developed by 
Ashby(Ashby & Cebon, 1993). 


The present study investigates the effect of filament orientation and AM-induced 
heterogeneous anisotropic behavior of 3D-printed structures on their performance under intrinsic 
shrinkage loading for 7 days determined by quantifying their physical characteristics (e.g., strength, 
fracture toughness) subject to quasi-static three-point flexural test, and compare them with of their 
cast counterparts. To assess how 3D-printed structure harness the heterogeneous interfaces to 
deflect the crack penetration through them, leading to damage localization at the restrained layer, 
X-ray microtomography (X-ray wCT) technique (1.e., no need of synchrotron facility) was used at 
two magnifications to analyze the microstructural characteristics of intact (1.e., not tested) damage- 
free and shrinkage-induced damaged 3d-printed and cast structures. In parallel, the digital image 


correlation (DIC) technique was employed to speculate the data from X-ray wCT on the evolution 
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of fracture zone size and microcrack opening associated with the shrinkage microcracking. 
Another goal of this study is to discern whether bioinspired Bouligand structural design as the 
second layer of resistance to bending can contribute to promoting a distinct mechanism by 
delocalizing the microcrack evolution in 3D-printed brittle hcp structures. A numerical framework 
is developed to assess the contribution of carefully designed structural orientations (found in 
biological systems) in offering a second layer of protection, thus enhancing shrinkage cracking 


resistance of the hcp element. 


5.3 Materials and methods 
5.3.1 3D Printing Setup 


Filament-based direct ink wiring (DIW) allows for the fabrication of ceramics utilized with 
a prudently designed hierarchal structure at multiple length scales based on extruding the 
pressurized ink without any heat application(Lewis, 2006). The filament-based approach to direct 
ink writing of 3D cement paste structure is based on extrusion freeform fabrication (EFF), where 
the ink is continuously extruded through a fine cylindrical nozzle to fabricate the filamentary 
structure. Figure 5-2 shows the integration of a stepper motor-driven ink delivery system 
(Structur3d Discov3ry Paste Extruder) and a gantry-based 3D printer (Ultimaker 2 Extended+) 
assembly to extrude the ink (henceforth, fresh cement paste) as a single continuous array of 
filament. The filament diameter (size) mainly depends on three printing factors of nozzle diameter, 
ink rheology, and printing speed. Given the favorable high yield stress and viscosity of printing 
paste, the DIW technique fabricates the structure of filaments that is capable of holding their 
designed shape. For further details on experimental configuration, the readers are referred to the 


conference paper by Moini et al. (Moini et al., 2019). 
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Figure 5-2. Filament-based extrusion system: (a) a stepper motor-driven extrusion system 
(Discov3ry) and a gantry-based 3D printer ; (b) constant-displacement extrusion through a 
syringe with lever pressure force at a uniform volumetric flow rate (adapted from Moini et 
al. (Moini et al., 2019)). 


5.3.2 Mixture Proportions, Mixing Procedure, and Curing Conditions 


For two purposes of structural integrity, including printability (1.e., extrudability) and 
buildability G.e., shape-holding), an iterative ink design process is established to introduce two 
favorable fresh properties of viscosity and yield stress. To this end, a sequential combination of 
viscosity modifying admixture (VMA) and superplasticizer additives are utilized to increase the 
resistance to segregation and facilitating the flowability of ink of 0.275 w/c, associated with a solid 
fraction of 53%. The final ink is composed of commercially available Type-I Portland cement 
powder, deionized water, HRWRA, and VMA. For every 250.0 grams of cement, 65.3, 1.84, and 
2.38 grams of deionized water, HRWRA and VMA are used, respectively. 


A Twister Evolution Venturi vacuum mixer is utilized to achieve a homogenous ink of air 


void-free, thus inhibiting the degradation of the hcp structural element. After a sequential 
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combination of water with HRWRA and VMA, the fluid-type mixture is added on the cement 
powder, and then the cement paste is mixed at a pre-mixing mode (of slow speed) for the duration 
of 25 s. Afterward, a torque of 400 rpm at a 70% vacuum level is applied to the paste for a duration 
of 90 s, which is followed by a final mixing process of applied torque of 400 rpm at a 100 % 
vacuum level for the duration of 90 s. Upon mixing, the final ink is loaded into the syringe, which 
is equipped with the plunger and assembled on the stepper-motor extruder. 

Extrusion of the 3D-printed specimen and placing cast specimens in the forms are carried 
out in a lab environment at 20 + 3 °C and 60 + 3 % RH. To apply restraint to these specimens 
while hindering any possible moisture exchange between the fresh specimen and the substrate (on 
which specimens are placed), we have utilized a novel approach of sanding Teflon rectangular 
prism of 50 mm thickness with a sandpaper grit of 80 to create a superhydrophobic surface 
possessing an RMS roughness of 14.5 um and a high advancing contact angle of 146° (Nilsson et 
al., 2010). Sanding is performed in an arbitrary manner for 20 s in any specific direction, followed 
by cleaning the sanded surface with acetone, rinsing the residues with deionized water, and 
cleaning with pressurized air. Figure 5-3 outlines the scanning electron microscopy (SEM) scans 
from untreated Teflon surface and Teflon surface sanded with 80-grit sandpaper, highlighting an 
introduction of surface roughness in sanded Teflon surface. Conversely, a thin Teflon sheet of a 
thickness (of 2 mm thickness) adhered to an acrylic rectangular prism of 10 mm thickness was 
built to fabricate a perfect superhydrophobic substrate with no roughness (1.e., unrestraint 
substrate). The former method of substrate fabrication is used to exert shrinkage-associated local 
tensile stress at the specimen-substrate interface and induce shrinkage microcracking along with 
the depth of “damaged” specimen. In contrast, the latter method is used to obstruct any tensile 
stress development at the specimen-substrate interface which leads to making a perfect 
“undamaged” specimen. Upon the end of specimen fabrication, the specimens are placed in a 
curing chamber of 20 + 3 °C and 33 + 2 % RH (using a saturated solution of Magnesium Chloride 
(MgClz)). Cast specimens are demolded at the age of 24 h and placed back in the chamber for the 


remainder of the experiment. 
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Figure 5-3. SEM images of Teflon surfaces sanded with 80-grit sandpaper: (a) smooth at 100x, 
(b) smooth at 200x, (c) smooth at 1000x, (d) 80-grit at 100x, (e) 80-grit at 200x, (f) 80-grit at 
1000x. 


5.3.3 Slicing and Design Parameters 


Simplify3D is the commercially available slicer that has provided detailed information on 
converting digital 3D models to 3d-printable parts, required amount of ink, printing time, and 
creating optimized printing toolpath. All processed information is generated in g-code which 
includes: Point cloud coordinates (X, Y, Z axes) of the toolpath to govern the travels of the nozzle 
and the bed; Extrusion to govern the quantity of extrusion with regards to the movement of the 
nozzle; Printing speed to control the speed of the movement of the nozzle. Simplify3D controls 
the E and F axes using an extrusion rate multiplier and printing speed in the slicer interface. 
Additional motions of the print head for moving the nozzle to ‘home’ coordinates without running 
into the printed objects are also programmed into the g-code command. 

The mere dimensions of the structure of the hcp element of this study is a rectangular prism 


of 50 mm x 12 mm x 12 mm. To achieve the desired lamellar and Bouligand architectured elements, 
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the geometrical characteristics of filament height, filament width, and structural orientation at each 
layer, combined with the processing parameters of printing speed, extrusion multiplier, and infill 
percentage are programmed in the slicer. Each filament layer 1s composed a continuous printing 
toolpath with two structural orientation characteristics of 0 and y, as shown in Figure 5-4; 6 defines 
the orientation of printed filaments in the first layer concerning the X-axis, while y identifies the 
re-orientation of printed filaments concerning their neighboring bottom layer along Z-axis. Figure 
5-4 (al-a3) depicts the sequential printing toolpath of the hcp element with the lamellar 
architecture of 6 = 0° and y = 0°, where Figure 5-4 (a4) highlights the arrangement of filaments 
which are prone to cracking under three-point flexural loading along the X-Z plane. Figure 5-4 
(b1-b3) depicts the sequential printing toolpath of an hcp element with the lamellar architecture of 
0 = 90° and y = 0°, where Figure 5-4 (b4) highlights the arrangement of filaments which represents 
the susceptibility of the interface region between filaments to cracking under three-point flexural 
loading along the X-Z plane. Figure 5-4 (cl-c3) depicts the sequential printing toolpath of an hcp 
element with Bouligand architecture of 6 = 90° and y = 7°, where Figure 5-4 (b4) highlights the 
arrangement of filaments which represents the susceptibility of twisting interface region between 
filaments to cracking under three-point flexural loading along the X-Z plane. A layer height (.e., 
filament height) of 1.0 mm and the extrusion width of 1.63 mm are programmed. The speed of 


printing is 750 mm/min. 
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Figure 5-4. Schematics of printing path and structural orientation of (al-a4) lamellar 
architectured structure, (bl-b4) lamellar architectured structure, and (cl-c4) bio-inspired 
Bouligand architectured structure : (al,bl,cl) Printing path of the first layer of filament 
highlighting the individual filament orientation concerning primary printing direction (along X- 
axis), 0 = 0°; (a2,b2,c2) Printing path of the sequential second layer of filaments highlighting 
the individual filament reorientation respect to of first layer filaments; (a3,b3,c3) resulting 3D- 
printed structures; and (a4,b4,c4) arrangement of filaments in the lamellar and Bouligand 
architectured hcp elements showing potential cracking paths under mechanical loading. 





5.3.4 X-ray uCT Imaging 


Because the interior microcracking development associated with inherently shrinkage 
behavior of cementitious material is unsightly on the surface of the element, the level of stress is 
significantly dissimilar inside the element experiencing inelastic deformation than that on the 
surface. Therefore, it is vital to observe the damage gradient along with the depth of 3D-printed 
and cast hcp elements. To this end, we employ X-ray wCT (micro-computed tomography) imaging 
technique to incorporate the logging of sequences of 2D radiographs recorded at numerous 
viewpoints around a revolving element. The entire volume (3D appearance of the interior and 
exterior) of the element is digitally rendered via mathematical reconstruction of re-configured 
spatial maps. Subsequent 3D renderings are characteristically displayed as a sequence of sliced 
radiographs with intensities associated with X-ray absorption and specific mass of material at a 
certain voxel (Landis & Keane, 2010). The subsequent disparities in intensity enable 
documentation of different phases and features of the intact and damaged structure at micro-level 


(e.g., microcracking) and their associated 3D distributions. Traditionally, X-ray microscopes 
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utilized for the uCT characterization tool utilize a flat panel detector and then depend on single- 
step (1.e., geometric) intensification, which implies the degradation of imaging resolution as the 
size of the specimen and the working distance increase. 

In this study, an X-ray microscope (XRM), Zeiss Xradia 510 Versa, is used to characterize 
the microstructural feature of both undamaged and damaged elements, including bulk elements 
and microcracks, which permits for a rise in the resolving power of X-rays via dual-stage 
intensification procedure. Therefore, a favorable field of view (FOV) is determined to examine the 
middle volume of the specimen through the intensification procedure, which includes maintaining 
spaces between the source, detector, and element (similar to traditional wCT technique). Using a 
detector equipped with a scintillator and objective lens allows us to achieve the highest resolution 
of the region, which is hypothesized to contain a high density of microcracks, resulting in 


degradation of the overall strength of the element under 3-point flexural loading. 


5.3.5 Multi-cutting compliance tests 


To quantitatively measure the length and nature of the microcrack growth behind the crack 
tip, a multi-cutting compliance technique, similar to that of Wittmann and Hu(Hu & Wittmann, 
1990), was employed for undamaged and damaged 3D printed and cast elements. Particularly, a 
thin diamond saw blade was used to extend the cut length consecutively and therefore 
incrementally abolish the crack wake (steps of ~0.6 - 1 mm), while two specimen compliances 
obtained from slopes of load-line displacement vs. load curve and crack mouth opening 
displacement (CMOD) vs. load curve were quantified (via the digital image correlation cut width 
(~ 100 um) exerts an insignificant effect on the compliance. As part of the crack wake, which is 
abolished, is free of traction, no alteration in compliance is observed after cutting (1.e., the cut 
region is free of shrinkage-induced microcracks). Though, when effective microcrack bridges are 
abolished from the crack wake by the diamond blade, there is a rise in compliance. Through the 
employment of this approach, the length of shrinkage cracking advancement along the specimen 
height is assessed by logging the length of the notch when the specimen compliance with active 
shrinkage cracking occurs to raise until it becomes identical to that of undamaged specimen 


compliance. 
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5.3.5.1 Load-line displacement vs. load compliance measurement on notched beams 


A quasi-static three-point bending test is carried out by MTS Universal Testing Instrument 
to measure the applied displacement and corresponding load. Before the specimens were notched 
on the diamond saw, load-displacement tests are performed to calculate the “un-notched” 
compliance, Co. The loading rate is justly fast enough (a displacement-controlled loading rate of 2 
umls) to reduce the creep effect. 

The value of critical fracture toughness, Kr can be calculated from the dimensions of the 
beam, the peak load, and the notch depth. Several equations are obtained by different stress- 
analysis technique but all effectively equivalent, appear in the literature; this work uses the one 
deduced by ASTM E-1820 (Standard Test Method for Measurement of Fracture Toughness: 
ASTM E1820 - 06, n.d.) as shown in equations 5-1 and 5-2: 

_ | Fmax> 
ie= 


ar f(a/w) 5-1 


where: 


2(1+26)(1-f) 


f(a/w) = 5-2 


In equation 5-1, Pmax indicates the peak load, Sis the span length, B is the average thickness, 
Bn is the net thickness, A is the beam height. In equation 5-2, a is the notch depth. 
The change of compliance with notch depth can be derived by making use of a general relationship 
between the critical strain energy release rate, Gr, and the specimen compliance(Kenny & 


Campbell, 1968), shown in equation 5-3: 
1 
Gic = 2 (Prax /B)* (dC /da) 5-3 


The stress analysis in an hcp element is analyzed based on the linear elastic fracture 
mechanics (LEFM) approach. Under linear elastic behavior before catastrophic failure assumption, 
the energy flowing through the fracture process zone (FPZ) per new crack surface growth is 


Kj-(1-v?) 


equivalent to Griffith’s energy release rate. Given Gre = , equations 5-1 and 5-3 are 


combined to obtain a relationship between C and a. To make the relationship generally applicable, 
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the elastic modulus must be eliminated from the equation; this 1s accomplished by calculating the 


ratio Co/C rather than C. Eliminating Kr and Gr and considering B = By leads to equation 5-4: 


dC  2(1—v*)S? 


— = —____ - 5-4 
== [f(a/h)] 
Hence 
2(1—v*)S* f [f(a/h)]? 
= — | = 5-5 
C T f n3 da 
Now, from elementary beam theory, 
ET 
E = ——— 5-6 
Co 4bh? 


where Cois the compliance of the beam when a = 0. The substitution of equation 5-6 into equation 

5-5 leads to equation 5-7: 
C a/h)]* 
oo 4 { DÈ 
C 

Where A = 


5-7 
8(1-v2)s7h Sia 
——.—\. For the geometry of the present experiment, specimen, A = 1.296 if v = 0.2. 


2 2 10 
Now, prem’ da = É (=) +e = (=) + constant = f(a/h) + constant. Since C/Co = 


1 when a = 0, C/Co =Af(a/h) + 1. 


5.3.5.2 CMOD vs. load compliance measurement 


The compliance method is a primary technique to determine the length of propagated crack; 
nevertheless, it is known that the presence of microcracks ahead of clearly cut notch influences the 
accuracy of the method significantly. Therefore, the effect of microcrack zone extension on overall 
compliance of specimens should be examined using a compliance calibration curve. This curve is 
constructed with cut specimens within their linear-elastic region. It is sensible to consider that the 
calibration curve is constructed for undamaged specimens, where no micro-damage associated 
with shrinkage ahead of clearly cut notch exists. Figure 5-5 shows a qualitative representation of 
measurements of loading compliance of the undamaged and damaged specimens (Cz and Cv) to 
examine the existence of crack bridging along with the clear-cut notch. Figure 5-5(a) outlines the 
experimental configuration of a three-point bending test applied to cement paste specimens of 


undamaged and damaged associated with shrinkage. A diamond saw is used to consecutively cut 
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out the specimen (and potential microdamage ahead of notch tip) at a specific increment to 
quantitively determine the extent of the microcracking zone behind the notch tip, given that a 
damaged material can be considered as a component constituted of materials of two elastic moduli 
of undamaged material (£) and material with microcracks (£*). Consequently, there is a difference 
in compliance due to the elastic modulus reduction shown in Figure 5-5(b) where P is the load, 
and CMOD is the crack mouth opening displacement. Figure 5-5(c) depicts the compliance 
calibration curve that is constructed with a notched specimen within their linear-elastic region. 
This figure illustrates how the length of microcrack bridging can be determined with the 
compliance method in conjunction with multi-cutting, where C specifies the compliance 
calibration curve (solid line), and Cp represents the compliance of a damaged specimen acquired 
from the multi-cutting tests. The x-axis in Figure 5-5(c) shows the increase of notch length with 
the origin at the initial notch length of ~0.6 mm. With the incremental cutting, the existence of 
bridging crack can be separated from the bridge-free region. When the part of the crack wake that 
is abolished contains bridging cracks, the compliance will be increased by cutting. When the 
cutting is through the entire crack wake region, the compliance of damage specimen and 
calibration curve will merge. as represents the initial notch size in this study, and aris the total 


propagation, including the extension of the active bridge zone. 
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Figure 5-5. Schematics of multi-cutting compliance experiment (a) investigation the effect 
of notch size extension on the mechanical response of undamaged and damaged cement 
paste elements under three-point bending test, (b) a loading and subsequent unloading 
response within linear-elastic region where microcrack zone length to be quantified, (c) 
evaluation of compliance curves from multi cutting of undamaged and damaged 
specimens. 
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5.4 Results 


As noted above, the interfaces between filaments in biological materials with the hierarchal 
pattern are one of the important structural features that determine the anisotropic behavior 
governed by the preferred alignment of filaments. Here, we emphasize the key role of interfaces 
in deflection and arrest of microcrack development associated with restrained shrinkage, which 
contributes to the promotion of a damage localization mechanism. This hypothesis is derived based 
on the fact that the interfaces are characteristically weaker than filaments in biological materials 
(W. Yang et al., 2019). In other words, the easy crack growth across the interfacial region 1s serious 
for rendering the toughening effect via the complexity of filament orientation. Besides, a carefully 
designed structural orientation of filaments plays a key role in allowing for developing an 
additional number of delocalized microcracks distributed throughout the domain subject to 


external loading before the catastrophic failure of the damaged hcp element. 


5.4.1 Mechanical performance of 3D printed hcp elements with lamellar architecture and 
cast hcp elements 


Here we examine the effect of shrinkage-induced internal micro-cracks on the mechanical 
performance of 3D printed hcp elements with two structural orientations. To characterize the 
strength of two main constituents of filament and interface in 3D printed elements, the values of 
modulus of rupture (MOR) of printed elements with two structural orientations of 0 = 0° and y = 
0° and 0 = 90° and y = 0° are quantified under three-point bending test. A comparison of the mean 
value of five MOR measurements for damaged 3D printed specimen (1.e., 3DP (D)) and damage- 
free 3D printed specimen (1.e., 3DP (UD)) highlights a reduction of 11.6%. In contrast, the value 
of MOR for damaged cast specimen (1.e., cast (D)) is reduced by 30.6% when compared to that 
for damage-free cast specimen (1.e., cast (UD)). The higher decrease of MOR in cast specimen 
leads to the idea that microcracks associated with shrinkage are easily distributed through its 
mesoscopic homogenous structure along with the specimen depth, effectively compromising the 
resistance of the material to damage. On the other hand, the lower decrease of MOR in 3DP 
specimens promotes a novel damage mechanism where the heterogeneous interface region 
contributes to trapping the microcrack development at the restrained filament layer. Owing to the 
weakness of interface by 25% when compared to the filament, promotes the hypothetical idea that 


these interfaces can deflect and channel the propagating microcracks from the restrained layer 
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along these predefined paths. To this end, we surmise that the fabrication of weaker interfaces 
increases the ability of 3DP elements in guiding and deflecting cracks to introduce a novel 
toughening mechanism through damage localization. The motivation of this work is inspired by 
biological systems where nature leverages from strong and brittle filaments delimited by weaker 
interfaces architectured in careful designs to induce the crack propagation along with the well- 


designed interface, effectively promoting damage and flaw resistance properties. 
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Figure 5-6. Mechanical response of damaged (D) and damage-free (UD) 3D printed (3DP) 
and cast hcp elements tested in three-point bending. 


5.4.2 Multi-cutting fracture toughness and compliance experiments 
5.4.2.1 5.4.2.1 Fracture toughness measurement 


The toughness represents the materials’ resistance to cracking, particularly by resisting the 
increase of clearly cut specimens. A prime source of reduction in fracture toughness is the 


development of inherent cracks formed at hierarchical length scales and distributed behind a crack 
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tip in brittle hcp material. When there are bridged or unbridged microcracks distributed at or behind 
the crack tip in damaged hcp elements (3DP or cast), the elements at the crack wake require lower 
energy to fracture when compared to their undamaged counterparts. Here, we present experimental 
results that a more quantitative description of how far shrinkage-induced microcracks propagate 
from the interface between specimen and substrate (1.e., the bottom of the specimen). Figure 5-7 
outlines the Kre: measurements of two 3D printed specimens of 3DP (UD) and 3DP (D) along with 
two cast specimens of the cast (UD) and cast (D) concerning various notch length, demonstrated 
by a/h ratio. The initial a/h ratio is 0.05 (0.05 x 12 = 0.6 mm), highlighted as Zone I. A comparison 
of findings from cast specimens indicates the presence of region with significant weakness behind 
the crack tip that is impacted by shrinkage microcracks. It is observed a reduction in Kr by 33% 
between a/h ratio of 0.05 and 0.163, whereas no significant deviation is observed at higher a/h 
ratios. Therefore, it is clear that the length of crack growth in cast specimens is extended to 2 mm 
(= 0.163 x 12 mm) during 7 days of exposure to the dry environment. A similar comparison of the 
Kı trend in 3D printed specimens highlights the limitation of the length of crack growth to 0.6 mm 
(= 0.05 x 12 mm), where a considerable deviation of 16.7% is observed. This implies that the crack 
propagation tends to propagate along with weak interfaces, requiring little energy, and inhibited to 
penetrate the neighboring top layer of filaments. The carefully designed interfaces are introduced 
within the element to alter its mechanical response to inherent shrinkage behavior entirely. On 
propagation a/h increases (> 0.05), no significant deviation of Kris observed due to the ability of 
(porous) weak interfaces for stabilizing cracks, effectively alleviating the strength compromise of 


the 3D-printed elements when compared to their cast counterparts. 
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Figure 5-7. Kr evolution of damaged and damage-free 3D-printed and cast hcp elements tested 
under three-point bending in the multi-cutting experiment. 


5.4.2.2 Compliance measurement 


The relationship between notch depth and load-line compliance, defined as beam 
displacement at the loading point divided by the load, of multi-cutting experiments for four 
specimen types of the cast (UD), cast (D), 3DP (UD), and 3DP (D) are shown in Figure 5-8. The 
theoretical relationship between Co/C and the notch depth is superimposed on the experimental 
results. The comparison of experimental results of each specimen type with the theoretical line 
indicates an insignificant deviation between the theoretical line and damage-free specimens (cast 
(UD) and 3DP(UD)), which implies that the theoretical line adequately represents the change in 
compliance with notch depth for all shallow notches. Conversely, the experimentally measured 
compliance values for damaged specimens (cast (D) and 3DP (D)) are lower than that for the ideal, 
undamaged specimen, providing direct evidence to the notion of microcrack extension at and 
behind the notch tip. Furthermore, Figure 5-8 shows that an essentially significant increase in the 
cast (D) and 3DP (D) is observed at the notch length of smaller than ~0.6 mm, highlighted in dark 
red color (Zone I). As the saw-cut notch length (1.e., microcrack wake region) is reached within 
~2mm (Zone II), highlighted in light red color, no increase in 3DP (D) specimen compliance is 
observed. In contrast, a significant variation is shown for cast (D) specimen compliance. This 


indicates that the microcrack extension length is on the order of 0.6 mm for 3DP (D) specimens, 
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while a greater value of 2 mm is found for cast (D) specimens. This experimental finding of load- 
line vs. load compliance obtained from the multi-cutting approach is in a strong agreement with 
that of fracture toughness measurement as they promote the effectiveness of interfaces in trapping 
the microcrack growth associated with shrinkage, and channeling the microcracks into stable well- 
designed configurations in the first “restrained” layer of filaments. The cast specimens are, 
however, deprived of this structural feature for controlling crack growth, which leads to being 


more likely susceptible to undergo excessive damage associated with shrinkage. 
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Figure 5-8. Relationship between compliance ratio and notch depth. C = compliance, Co = 
compliance for zero notch. 


As has been noted in this work, the present results affirm the notion that the shrinkage- 
induced crack localization is a major competitive mechanism in 3D-printed specimens. New 
results on mechanical performance characteristics of 3d-printed and cast specimens confirm how 
3D configuration of interface regions provide ideal crack propagation sites, effectively inhibiting 


crack growth into the neighboring layers. The load-line vs. load compliance measurement using 
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multi-cutting compliance experiments showed how the greater extension of microcracks along 
with the specimen height leads to a significant increase in the global compliance of the damaged 
specimen when compared to its undamaged counterpart. Here, we perform a compliance-based 
measurement of crack mouth opening displacement (CMOD) of clearly saw-cut notch regarding 
the applied load to analyze the effect of microcrack length extension on local compliance of 
specimen. To determine the CMOD in a non-contact optical mode, digital image correlation (DIC) 
is used to extract the real-time full-contour displacement on a specimen surface(Yates et al., 2010). 
For more information on specimen preparation and image recording method, readers are referred 
to an article by Sutton et al. (Sutton et al., 2007) examining the capability of the DIC method for 
quantification of deformation and CMOD in ductile aluminum subjected to mixed-mode loading. 

Representative CMOD-load responses are shown in Figure 5-9 for the undamaged and 
damaged cast specimens, which demonstrated the local compliance alteration about sequential 
notch length extension ranging from a = ~0.6 mm to ~4 mm. The line fitted to the data points 
suggests a best-fit linear function to determine the compliance value approximately. Shown for 
comparison at a specific notch length, it is apparent that the experimentally acquired compliance 
values of damaged specimens (in Figure 5-9(b)) are greater than those for the ideal, undamaged 
specimens (in Figure 5-9(a)) at notch sizes as long as roughly 2 mm. When further extending the 
notch length size to larger values (1.e., 3 mm and 4 mm), the difference of compliance values of 
damaged and undamaged cast specimens appears to be negligible, which implies a microcrack- 


free zone at these notch length sizes. 
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Figure 5-9. Quantitative representation of pre-failure CMOD-Load curves for (a) 
undamaged cast specimens with notch length range between a = 0.6 mm and a = 4 mm, 
and (b) damaged cast specimens with notch length range between a = 0.6 mm and a = 4 
mm. 
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Similar multi-cutting experiments and analysis are performed on undamaged and damaged 
3DP specimens, shown in Figure 5-10. A comparison of compliance values of undamaged and 
damaged specimens at a specific notch length indicates that microcrack extension zone for 3DP 
specimen is significantly smaller than that for cast specimen; indeed, the zone length is on the 
order of 0.6 mm long in damaged 3DP specimens when compared to their undamaged counterparts. 
It is concluded that the microcrack zone is extended to a length on the order of ~ 2 mm, which is 
much larger than those inferred for cast specimens that yielded a steady state “plateau” compliance 
after ~ 0.6 mm. This difference is seemingly due to the fabrication of weaker interfaces within the 
bulk material to induce the crack propagation along these interfacial regions, thereby providing a 


layer of protection against detrimental shrinkage cracking. 
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Figure 5-10. Quantitative representation of pre-failure CMOD-Load curves for (a) 
undamaged 3DP specimens with notch length range between a = 0.6 mm and a = 4 mm, 
and (b) damaged 3DP specimens with notch length range between a = 0.6 mm and a = 4 
mm. 
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Figure 5-10 continued 
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The efficient microcrack advancement determined from the compliance-calibration curve 
is used. The associated compliance curve for both cast and 3DP specimens obtained from multi 
cutting and highlighted with black color is presented in Figure 5-11(a) and (b), respectively. Each 
data point represents the mean value of five measurements for the associated notch size, while the 
curve fitted to the data points specifies a best-fit sixth-order polynomial function to the data. It is 
seen that the variation of the compliance through the notch length increase in damaged specimens 
is less than that of the compliance-calibration curve (1.e., compliance variation of the undamaged 
specimen) since the bridging stress transported within the zone with shrinkage-induced 
microcracks is less than what is likely to be transported within the zone with intact material. 
Through the curve fitting, it is discovered that the lengths of the microcrack zone extended along 
the specimen height are on the orders of 2.57 mm and 1.23 mm, which is in agreeable compliance 


with previous findings. 


154 


0.1 
0.09 





A Z 0.08 
S S 
£ = 0.07 
= 2 0.06 
as ae 
S L 0.05 
5 5 
Z = 0.04 
E- & 0.03 
5 5 + 18.7% + 5.8% rA” 
O © 0.02 @-5-9) 
© 
ui \ 
3DP (UD) 
0 
0 | 2 3 4 
Notch length, a (mm) Notch length, a (mm) 
(a) (b) 


Figure 5-11. Results from multi-cutting compliance experiments on (a) an undamaged and a 
damaged cast specimen, and (b) an undamaged and a damaged 3D-printed specimen. The 
indicated location of the crack tip is measured by multi-cutting methods, while each plotted data 
point corresponds to the mean of five compliance measurements. 


5.4.3 Micro-imaging of crack growth 


This work investigates the effect of interface characteristics on alleviating the growth of 
shrinkage-induced damage within the cast and 3D printed lamellar elements using the X-ray wCT 
scanning method. Traditionally, laboratory-based X-ray radiographs (1.e., images) are produced 
from an X-ray source bombarding high-speed electrons to the target region of interest (ROI) of an 
element (Landis & Keane, 2010). The resulting alteration in distributed intensities depending on 
the X-ray absorption and material density allows for the identification of various phases and their 
geometrical properties. 

In this study, an X-ray microscope (XRM), Zeiss Xradia 510 Versa, is used, which allows 
for scanning resolution enhancement through a dual-stage magnification technique. Figure 5-12 
depicts an overview of the X-ray wCT scanning procedure. Given the internal crack opening size 
associated with the shrinkage is on the order of microns, the microstructure of elements is 
characterized using 4X magnification with a resolution size of ~15 um. The desirable region of 
interest (ROD) is selected to enable us to visualize the randomly distributed microcracks within a 
cylinder of 12 mm in diameter and 12 mm in height located in the center of elements. A beam 


energy of 150 KeV, a power of 10 W, exposure time of 4 seconds, and full 360° rotation are used 
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for 4X scans of elements. For the post-processing of the images, Dragonfly software is used. One 
intact cast, one intact 3D printed, one damaged cast, and one damaged 3D printed prismatic hcp 


elements (50 mm x 12 mm x 12 mm + | mm) are used in this experiment. 
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Figure 5-12. Schematic illustration of procedure of X-ray uCT scan: (a) a laboratory-based 
X-ray technique generating and reconstructing cross-sectional images, (b) experimental 
configuration of specimen scanning using an X-ray microscope (XRM), (c) 3d rendering 
of the volume of the region of interest in 4X scan mode. 


5.4.3.1 Restrained and unrestrained shrinkage in 3D printed specimens with lamellar 
architecture 


To investigate the characteristics of the interface and its role in trapping the microcrack 
growth, the microstructure of 3D printed hcp element with 0 = O° is characterized, and the 
microcrack distribution at carefully-selected cross-sections is studied, as shown in Figure 5-13- 
Figure 5-15. In wCT images of hcp, darker intensities represent both pores filled with air or water 
and microcrack openings, with greyscale intensities associated with hydrated cement paste 
particles. 

Three horizontal slices (F1, F11, Fiz) and three horizontal slices (lı, lio, 11) shown 1n cross- 
sectional view in Figure 5-13(a), corresponds to the ‘core’ (i.e., through the interface) and 
‘interfacial regions’ between the filaments. It 1s necessary to note that the numbering starts from 
the bottom of the specimen. Moreover, the horizontal slice (S) and two vertical slices (Hi, H2) 
corresponds to the ‘exterior’ cross-sectional views of the specimen that are exposed to exterior 


boundary conditions (1.e., slice S corresponds to the bottom surface of the filament in the first layer 
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deposited on the substrate with no roughness, whereas slice Hı and H2 represents the stack of 12 
layered filaments exposed to severe low humidity). Analysis of Figure 5-13(b1-b7) indicates no 
significant microcrack onset or growth within the filaments at the bottom layer. In contrast, 
observable microcrack distribution is observed within the filaments located at the two top layers 
(Fi; and Fi2 cross-sectional views) and the interface between them (Ii; cross-sectional view) owing 
to the drying shrinkage mechanism exacerbated by the severe dry ambient environment. 
Furthermore, observable microcrack growth in an arbitrary pattern on two lateral sides of the 
element (displayed in Figure 5-13(b8-b9)) provides further direct evidence of the deleterious effect 
of drying shrinkage on generating damage within filament microstructure. While a previous study 
has suggested for appropriate curing for alleviating the excessive degree of shrinkage cracking in 
3D printed elements, the comparison of MOR of two unrestrained cast and 3D printed hcp elements 
indicates ‘zero’ contribution of these developed micro-cracks on the mechanical performance of 


3D printed elements. 


157 


A 
A 


umi 768°6 
tal cep TT 


á 


~] 
Yə 
Qo 
~— 
— 
5 


a 


T=A 


um [¢0'T 
wr 7L8‘01 


~ 
on 
ey) 
neat 
SS 
ON 





Figure 5-13. X-ray uCT images of the microstructure of 3D-printed solid prism (0 = 0°) 
specimen containing shrinkage-induced microcracks collected during 4X scan: a) 2D image 
of the cross-section in XY plane, b1-b10) 2D images of the cross-section, highlighted in (a), 
in XZ plane, screening the shrinkage-induced microcrack propagation along with the 


specimen height. 


Here, the cracking behavior and associated damage mechanism in restrained 3D printed 
hcp elements are investigated to highlight the role of the interface in governing the damage 
mechanism using wCT images at 4X magnification. Figure 5-14 depicts the microstructural 
characteristics of the element at nine critical cross-sectional views in X-Z plane (slice S 
corresponds to the bottom surface of the filament in the first layer deposited on the substrate with 
a specific roughness; five slices Fi, F2, F3, F11, Fi2 corresponds to the ‘core’ (1.e., through the 
interface); and, four slices Ij, I2, lio, Iii are associated with the ‘interfacial regions’ between the 
filaments). These slices are intentionally selected to monitor the guiding capability of interfaces 
along their predefined surface (1.e., path). When the five bottom horizontal cross-sectional views 
are examined, a gradient of microcrack population and distribution is observed along the filaments 
of the first layer, with a maximum degree of cracking screened on slice S (1.e., Y = 0 mm) and a 
minimum degree of cracking screened on slice I (1.e., Y = 839 mm). When the second layer of 


filaments are examined microscopely, ‘zero’ cracking population and distribution is observed 
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which reveals the channeling of cracks into weaker heterogeneous interfaces and inhibiting the 
damaging of stronger filaments at the neighboring upper layers (slices F2, F3) and the interface 
between them (slice Io). Similar to unrestrained 3D printed hcp element, the microcrack 
advancement within filaments of the top two layers (slices F11, li, F12) is observed owing to the 
excessive degree of drying shrinkage. Following the influence of these pronounced interfaces on 
deflecting the overall crack paths, no crack population or distribution been observed on cross- 


sectional lı view, highlighting the microcrack channeling along with the interface between two 


filaments. 
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Figure 5-14. X-ray uCT images of the microstructure of 3D-printed solid prism (0 = 0°) 
specimen containing shrinkage-induced microcracks collected during 4X scan: a) 2D image 
of the cross-section in XY plane, b1-b10) 2D images of the cross-section, highlighted in (a), 
in XZ plane, screening the shrinkage-induced microcrack propagation along with the 


specimen height. 


Furthermore, wCT characterization of restrained 3D printed hcp element along the Y-Z 
plane (1.e., third plane) is necessary to monitor the realistic three-dimensionality of the microcrack 
growth gradient along with the specimen height and associated damage mechanism. The uwCT 
analysis 1s performed as a mean to improve further our understanding the central role of the 
interface, as shown in Figure 5-15. Figure 5-15 depicts the microstructural characteristics of the 


element at nine critical vertical cross-sectional views in Y-Z plane (two slices Hı and H2 
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corresponds to the cross-sectional views of lateral filament layers directly exposed to the dry 
ambient environment, and seven slices H3-Ho correspond to ‘interior’ filament layers). Similar to 
unrestrained 3D printed hcp element, a significant degree of microcracking distribution is observed 
along the specimen height shown in slices H; and H2 (Figure 5-15(b1-b2)), where its onset occurs 
within the filaments located at the top layer propagating toward to the bottom layers. Figure 
5-15(b3-b9) depicts the examination of microcrack propagation within the stacked filaments in the 
Y-Z plane. The observational microcracking on each slice reveals the development of microcracks, 
initiating from the (restrained) bottom of element extended towards the interface between the 
filaments of the first and second layers with an associated maximum height of 919 um. According 
to wCT images in Figure 5-14, this microcracking zone length is approximately identical to the 
height of the first filament layer (~939 um), highlighting the influence of weak interfaces in 
hindering the crack propagation into the second layer of filaments. This finding complements our 
earlier observations that the mechanical response analysis that local and global compliances of 
damaged 3D printed hcp elements cut with a notch length up to 1 mm are higher than those of for 
undamaged 3D printed hcp elements. In other words, the fabrication of stiff and hard filaments 
delimited by weaker interfaces arranged in specific architectures can promote a unique damage 
mechanism that enhances the mechanical response of hcp elements without sacrificing their overall 


strength. 
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Figure 5-15. X-ray uCT images of the microstructure of 3D-printed solid prism (0 = 0°) 
specimen containing shrinkage-induced microcracks collected during 4X scan: a) 2D image 
of the cross-section in XY plane, b1-b10) 2D images of the cross-section, highlighted in (a), 
in XZ plane, screening the shrinkage-induced microcrack propagation along with the 
specimen height. 


5.4.3.2 Restrained and unrestrained shrinkage in cast hcp elements 


To highlight the importance of interface fabrication in enhancing the flaw-tolerance 
properties of 3D printed hcp element, the microstructure of its cast counterpart is characterized, 
and the microcrack distribution is studied at favorable cross-sections, as shown in Figure 5-16- 
Figure 5-18. 

Six horizontal slices (S, Li, L2, L3, L4, Ls) shown in cross-sectional view in Figure 5-16 
(a), corresponds to the mapping of microcracking population and distribution evolution along with 
the specimen height, while the two vertical slice (Hı, H2) corresponds to the ‘exterior’ cross- 
sectional views of the specimen that are exposed to the external boundary conditions. Analysis of 
Figure 5-16 (b1-b6) indicates the only significant microcrack pattern occurs at the top surface of 
the specimen (1.e., Y = 11,740 um), where the lateral and top surfaces of the element is subject to 
an excessive degree of drying shrinkage. When compared to its 3D printed counterpart, it is 
observed that less microcracking takes place in cast hcp elements deposited on unrestrained 
substrate since the formwork essentially plays a central role in decreasing the degree of drying 


shrinkage by limiting the exposure of the element to top surface. Similar to its 3D printed element, 
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observable microcrack growth in an arbitrary pattern is observed on two lateral surfaces of the 
element (displayed in Figure 5-16 (b7-b8)), which provides direct evidence of the deleterious 
effect of drying shrinkage on generating damage within filament microstructure. It should be noted 
that these microcracks have a less pronounced influence on the overall MOR of the hcp element. 

Six horizontal slices (S, Li, L2, Ls, L4, Ls) shown in cross-sectional view in Figure 5-16 (a), 
corresponds to the mapping of microcracking population and distribution evolution along with the 
specimen height, while the two vertical slice (H1, H2) corresponds to the ‘exterior’ cross-sectional 
views of the specimen that are exposed to the external boundary conditions. Analysis of Figure 
5-16(b1-b6) indicates the only significant microcrack pattern occurs at the top surface of the 
specimen (i.e., Y = 11,740 um), where the lateral and top surfaces of the element is subject to an 
excessive degree of drying shrinkage. When compared to its 3D-printed counterpart, it is observed 
that less microcracking takes place in cast hcp elements deposited on unrestrained substrate since 
the formwork essentially plays a central role in decreasing the degree of drying shrinkage by 
limiting the exposure of the element to top surface. Similar to its 3D-printed element, observable 
microcrack growth in an arbitrary pattern is observed on two lateral surfaces of the element 
(displayed in Figure 5-16(b7-b8)), which provides direct evidence of the deleterious effect of 
drying shrinkage on generating damage within filament microstructure. It should be noted that 


these microcracks have a less pronounced influence on the overall MOR of the hcp element. 
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Figure 5-16. X-ray wCT images of the microstructure of unrestrained cast prism specimen 
collected during the 4X scan: a) 2D image of the cross-section in the XY plane, b1-b8) 2D 
images of the cross-section along the XZ plane. 


Next, the cracking behavior and associated damage mechanism in restrained cast hcp 
elements are investigated to understand the governing damage mechanism using uCT images at 
4X magnification. Figure 5-17 depicts the microstructural characteristics of the element at eight 
critical cross-sectional views in X-Z plane (slice S corresponds to the bottom surface of the 
filament in the first layer deposited on the substrate with a specific roughness; eight slices Li, Lz, 
L3, L4 represent the microstructural details of the element adjacent to the restrained bottom surface; 
and, three slices Ls, Le, L7 represent the microstructural details of the element adjacent to the to 
surface subject to the excessive degree of drying). These slices are intentionally selected to monitor 
how damage-tolerant the hcp element is when exposed to an extreme drying environment, guiding 
the capability of interfaces along their predefined surface (1.e., path). When the five bottom 
horizontal cross-sectional views are examined, a gradient of microcrack population and 
distribution is observed along with the specimen height, with the highest degree of microcracking 
screened on slices S and L; and L2 (1.e., Y < 762 mm). As the slice location moves away from the 
bottom of the specimen (1.e., slices L3 and L4), the degree of microcracking reduces significantly 


due to the less pronounced tensile stress development. Similar to unrestrained cast hcp element, 


163 


the microcrack advancement within the portion of the element located adjacent to the top surface 
(on three slices Ls, Le, L7) is observed owing to the excessive degree of drying shrinkage. In 
overall, a more pronounced degree of microcracking with regards to microcrack population and 
distribution along with a higher microcrack propagation depth is observed in restrained cast hcp 
elements when compared to the 3D printed counterparts. This microcracking pattern eventually 
leads to a significant reduction of resulting mechanical response (e.g., MOR) of cast element, 


which suitably promotes our earlier findings in section 5.4.1. 
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Figure 5-17. X-ray wCT images of the microstructure of restrained cast prism specimen 
collected during 4X scan: a) 2D image of the cross-section in XY plane, b1-b8) 2D images 
of the cross-section along the XZ plane. 


Furthermore, uCT characterization of restrained cast hcp element along the Y-Z plane (..e., 
third plane) is essential to monitor the third dimension microcrack pattern gradient along with the 
specimen height. The wCT analysis is performed as a means to improve further our understanding 
of how the hcp element with no structural design of interfaces delimiting the bulk material is 
susceptible to the inherent damage associated with shrinkage, as shown in Figure 5-18. Figure 5-18 
depicts the microstructural characteristics of the element at eleven critical vertical cross-sectional 


views in Y-Z plane (two slices Hı and Hi; corresponds to the cross-sectional views of filament 
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layers neighboring exterior surfaces of the element directly exposed to the dry ambient 
environment, and nine slices H2-Hio correspond to ‘interior’ filament layers). Similar to 
unrestrained cast hcp element, a significant degree of microcracking distribution is observed along 
the specimen height shown in slices Hı and Hi; (Figure 5-18 (bl, b11)). In converse to its 3D 
printed counterpart, the microcracking pattern is observed to be distributed and extended across 
the entire microstructure. When examining the microcrack propagation height along nine interior 
slices shown in Figure 5-18 (b2-b10), it is revealed that the development of microcracks, initiating 
from the (restrained) bottom of element extended towards the specimen height reaches to a 
maximum length 1,580 um. When compared to uCT images of its 3D printed counterpart, this 
microcracking zone length is approximately ~ 1.7 times larger in restrained cast hcp elements, 
which contributes to sacrificing the resulting MOR of the element. This finding complements our 
earlier observations with regards to the mechanical response analysis that local and global 
compliances of damaged 3D printed hcp elements cut with a notch length up to 2 mm are higher 
than those of for undamaged 3D printed hcp elements. In other words, the absence of weaker 
interface fabrication leads to worsening of overall mechanical performance and damage tolerance 


properties of the hcp element. 
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Figure 5-18. X-ray wCT images of the microstructure of restrained cast prism specimen 
collected during 4X scan: a) 2D image of the cross-section in XY plane, b1-b11) 2D images 
of the cross-section along the XZ plane. 
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5.4.4 Finite element modeling and analysis: Bio-inspired architecture effect 


In this section, the mechanical response of restrained 3D printed cement paste elements 
with specific architectures subject to shrinkage loading, along with the associated damage 
mechanisms, have been investigated by examining and comparing the behavior of these elements 
with their unrestrained counterparts under flexural loading. Here, we focus on the Bouligand 
architecture inspired by the light-weight appendages of mantis shrimp (Figure 5-19 (a)) (Weaver 
et al., 2012). The last segment of this appendage is the dactyl club (Figure 5-19 (b)), which is 
composed of three main regions of impact (1), periodic (II), and striated (III). The architecture of 
the periodic region is found to follow a helicoidal arrangement of stacked layers of unidirectional 
chitin fibrils (Figure 5-19 (c)). A schematic illustration of an idealized Bouligand structure is 
shown in Figure 5-19 (d-e), where a pitch angle (y) is used for structural characterization. It is 
noteworthy that this architecture produces many crack deflection pathways during loading (Moini 


et al., 2018; Suksangpanya et al., 2017). 
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Figure 5-19. (a) Photograph of Mantis Shrimp. (b) Schematic image of the transverse 
section of the dactyl club (highlighted in Figure 5-19(a)) where I, II, and HI display the 
impact, periodic, and striated regions, respectively. (c) SEM micrograph of a fractured 
surface of the periodic region. (d-e) Representative prism creation: the multilayer section 
highlighted in multiple colors is built from the orientation of each layer with a pitch angle 


y. 
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To characterize the mechanical response of 3D printed hcp elements, we present a 
combination of simulation (3D finite element model), additive manufacturing, and flexural testing 
of Bouligand architecture to expound the role of multiscale architectures in hcp elements’ 
enhanced shrinkage-induced crack resistance. This work addresses the modeling of the behavior 
of 3D printed hcp elements through a two-step multiscale approach that first starts at the lower 
length scale, where the configuration of four semi-rectangular filaments. The second step considers 
the modeling of the entire 3-D printed prism hcp element composed of consecutive layers of 
filaments stacked on one another. Moini et al. (Moini et al., 2019) reported the self-drifting of 
filaments from their programmed toolpath, which leads to the creation of both diamond- and 
triangular-shaped gaps between neighboring filaments. In the first step, this work develops a three- 
dimensional finite element model- (FEM-) based model that accounts for the influence of this 
characteristic of printed microstructure on the overall tensile mechanical response of interface, as 


shown in Figure 5-20. The goal of this step is to homogenize the fracture response for a free-of- 


int,num int,num 
and Gro 


gap interface element in terms of cohesive properties (Tn ) to significantly 
minimize the computational cost required to discretize the curves generated during the filament 
deposition. To test the validity of fracture behavior of interface in this new configuration, a 
comparison between the overall mechanical responses of (1) four-filament configuration with a 
diamond-shaped gap and known tensile mechanical properties of both filament and interface 
phases (Figure 5-20(a)), and (11) four-filament configuration free-of-gap with known tensile 


mechanical properties of filament and homogenized tensile mechanical properties of the interface 


(Figure 5-20(b)) are performed subject to a mode-I loading condition. 
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Figure 5-20. Numerical framework representing the first step of multiscale modeling of the 
mechanical response of (a) a four-filament configuration with semi-rectangular shape in the 


filament axis, knowing the tensile strength of filament and interface (T J and T M and 


fracture energy of filament and interface iG and G E “€XP) obtained from experiment; (b) 
a four-filament configuration with rectangular-shape in the filament axis, knowing the 


tensile strength and fracture energy of filament (T A and Gi") obtained from experiment 


and tensile strength and fracture energy of the new interface (T{"°"™™ and Gi?’™™ ) 
calculated from numerical modeling. 


The elastic and cohesive properties of two phases of filament, as well as an interface in the 
four-filament configuration with the diamond-shaped gap (labeled as ‘int, exp’) obtained from the 
earlier experimental analysis, are outlined in Table 5-1. Given these properties along with the 
structural characteristics of two configurations, a trial-and-error process is employed until 
desirable cohesive properties for the free-of-gap interface (labeled as ‘int, num’) are obtained to 
calculate the homogenized stress-displacement (Gn-on) along with loading axis. This process 
includes first calculating the resulting peak stress obtained from the mechanical response of four- 
filament with diamond-shaped configuration and then calculating a desirable value for òn by 
calibrating the resulting On-6n of four-filament free of gap configuration to match with of that for 
its counterpart containing diamond-shaped gap, shown in Figure 5-21. Then, the desirable values 
for cohesive properties of the free-of-gap interface are measured and outlined in Table 5-1. Note 
that the intrinsic cohesive properties of filament and interface in pure mode-I are identical to those 


in pure mode-II. 
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Table 5-1. Material parameters of filament and interfaces in two various configurations. 
Mechanical 


ance Unit Filament (fil) Interface (int, exp) Interface (int, num) 
E GPa 6000 - - 
V - 0.2 - - 
Gic MPa.mm 0.00816 0.003 0.00267 
L nax MPa 10.5 8.81 6.57 
Ôn mm 0.001556 0.000681 0.000812 
Lez mm 0.45 0.115 0.186 
8 
—Diamond-shaped gap 
7 | —Free of gap 
6 
5 
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Figure 5-21. Validation of cohesive properties of four-filament free of gap configuration with 
of those for four-filament configuration with diamond-shaped gaps in terms of averaged stress- 
displacement response. The insets represent the evolution of the cracking of interface at 
corresponding loading applied to the configuration. 


While experiments are useful to quantify the mechanical properties, including strength and 
fracture toughness of hcp elements, they provide no knowledge of the competing damage 
mechanism in elements of complex (micro-)structures. To this end, we employ a FEM-based 


numerical approach to look at how the cracking onset and propagation takes place in anisotropic 
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heterogeneous composites which cannot be explicitly described by theoretical or experimental 
approaches. First, to verify the validity of the numerical model, the mechanical responses of 
printed hcp prisms with lamellar architecture are quantified using numerical analysis and then 
contrasted with earlier experimental observations. Following the experimental protocol, the three- 
point flexural loading in simulated using three rigid rollers having frictionless contact between the 
roller and the hcp prism (40 mm x 12 mm x 12 mm). The elastic properties of continuum elements 
are obtained from performing separate hcp prisms with no interface elements tested under identical 
loading conditions and then compared with the experimental data of unrestrained hcp elements. 
As a result, the values of elastic modulus E£ and Poisson’s ratio v are determined to be roughly 
about ~ 6000 MPa and ~ 0.2, respectively. 

Here, we quantify the o;-d response of unrestrained hcp element with two different 
characteristic angles 0 = 0° and y = 0° where loading plane is perpendicular to the filament axis 
(1.e., testing the mechanical response of the filament layers explicitly) as well 0 = 90° and y = 0° 
where loading plane is parallel with the interface axis (1.e., testing the mechanical response of the 
interface layers clearly), and compare them with experimental observations, as shown in Figure 
5-22. Comparing the numerical solution and experimental observation, it is seen that the 
percentage difference for calculation of peak stress becomes less than 4%, which implies a 
promising agreement. The next step is to predict the mechanical response of restrained hcp element 
with characteristic angles 0 = 0° and y = 0°, where it accounts for the contribution of shrinkage- 
induced microcracks propagated across the filaments of the first layer and neighboring interfacial 
region. Given the mechanical properties of intact and free-of-damage filaments in the above layers 
(1.e., layers 2-12), a trial-and-error process is employed until desirable cohesive properties for the 
filaments in the first layer are obtained to match with the mechanical response obtained from 
experimental data. After several trials and errors, the most desirable values of E, Tmax, Gic, and ôn 
are determined to be 3000 MPa, 5.1 MPa, 0.004 MPa.mm, and 0.00081 mm, which yield to a 
percentage difference of less than ~ 2% for peak stress value. Correspondingly, the favorable 
values of E, Tmax, Gic, and ôn for the damaged interface in the first layer are determined to be 3000 


MPa, 3.38 MPa, 0.0013 MPa.mm, and 0.00081 mm. 
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Figure 5-22. Mechanical response of lamellar architecture using the three-point flexural 
test: a) flexural stress-displacement (o;-0) for unrestrained (undamaged) and restrained 
(damaged) printed prisms with characteristic angles of 0 = O and y = O, as well as 
unrestrained printed, prisms with characteristic angles of 0 = 90 and y = 0. The insets 
represent the evolution of cracking at corresponding loading applied to the hcp element. 


In bio-inspired fibrous Bouligand architectures, a competing damage mechanism between 
solid failure, representing fiber fracture, and interfacial failure, representing the separation 
between fibers, is expected to control the trade-off between large and small pitch angles 
(Suksangpanya et al., 2017; Weaver et al., 2012). While small pitch angle allows the crack to 
deflect into the interface, the activation of this interfacial failure can be achieved by designing a 
crack path that facilitates the crack twisting. This provides delocalization damage mechanisms and 
the spread of damage to the neighboring interfacial regions. Therefore, this study focuses on the 
examination of the hcp elements with small pitch angles (e.g., y = 1°, 3% and 5°). This section aims 
to assess the use of Bouligand architecture as the second layer of protection for hcp elements when 
subject to damage associated with shrinkage. Three-point flexural testing is performed on these 
elements for the determination of resulting or-ò behavior and specific strength (MOR). 

The results presented in Figure 5-23 (a) illustrate that peak stress values consistently 


increases with the increase of pitch angle in unrestrained hcp elements (UD). At the same time, 
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independent of the architecture and structural orientation of filaments, a catastrophic failure of 
specimen takes place after reaching to the peak stress. Specifically, the observed increases in the 
peak stress are ~ 11% and 21% for architectures with y = 3° and y = 5°, respectively, when 
compared to their counterpart with y = 1°. While there is no available data on the effect of 
restrained shrinkage cracking on the mechanical performance of hcp elements with Bouligand 
architecture, we utilize the numerical tool to predict the resulting mechanical behavior and an 
integrated localization-delocalization damage mechanism through trapping the microcracking 
associated with shrinkage at the first layer of filaments, followed by guiding the cracks induced by 
mechanical loading through interfaces with stable configurations. When a comparison of 
6-0 responses of Bouligand architecture specimen with y = 3° obtained from numerical analysis 
and experimental data, a similar peak stress value is observed with a percentage difference of less 
than ~ 6%. Afterward, the cohesive properties of damaged filaments and interfaces produced in 
the first layer and interfaces between the filaments in the first layer and consecutive top layer (1.e., 
second layer) obtained from earlier numerical results are incorporated into the first layer of printed 
Bouligand architecture to mimic the microcracking associated with shrinkage. A similar 
6,-0 response between unrestrained (UD) and restrained (D) Bouligand architecture specimen with 
y = 3° 1s observed, except an insignificant increase in the specimen compliance and an insignificant 
decrease in the peak stress value are depicted. 

In terms of strength, the effect of structural anisotropy fabricated by additive manufacturing 
on alleviating the effect of restrained shrinkage cracking is quantified in three distinct hcp elements 
of the cast, printed lamellar architecture, and printed Bouligand architecture, as shown in Figure 
5-23 (b). The decrease in the MOR values of 3D printed element with 0 = 0° and y = 0°, 3D printed 
element with 0 = 90° and y = 0°, 3D printed element with 0 = 90° and y = 3°, and cast element are 
11.66%, 9.8%, 5.76% and 29.51%, respectively, highlighting the important influence of structural 
orientation of filaments on solving a primary durability issue. Furthermore, it is observed that the 
unrestrained printed Bouligand architecture with y = 5° and unrestrained printed lamellar 
architecture hcp elements are statistically identical, emphasizing on the role of crack deflection at 
the interface as the primary source of enhanced damage and flaw tolerance. 

The detailed images of typical crack propagation in a stepwise pattern, crack reorientation, 
and branching offering enhanced fracture properties for unrestrained and restrained printed hcp 


elements with y = 3° are illustrated at various loading stages in Figure 5-23 (c—n): 1) At the pre- 
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peak response before peak stress (highlighted in green color), the numerical modeling suggests 
that a particular primarily unidirectional crack propagates along with the vertical interfaces of first 
two layers in the Y-Z plane (shown in Figure 5-23 (c)). In Figure 5-23 (d) and (e), the crack path 
is detected to be redirected to the weak interface along the filaments (1.e., it is parallel to X- 
direction). This stepwise crack pattern promotes the idea of using Bouligand architectures in which 
they tend to allow the cracks propagate in twisted designs following the path of the interfaces. 
These twisting patterns have been explored to be responsible for increasing damage resistance and 
promoting the spread of damage (1.e., delocalization damage mechanism). As loading progresses 
passing the peak stress (highlighted in yellow color), the crack propagation is systematically 
deflected into interface channels rather than deflecting straight into the filaments, where maximum 
concurrent tensile and principal shear stresses are most likely to be exerted (shown in Figure 
5-23 (f), (g), and (h)). At certain stages of loading before ultimate failure (highlighted in purple 
color), it is observed that these weak interfaces are exploited to deflect and branch the cracks into 
the design where growth becomes more energy-consuming, which is the strategy observed in 
biological systems (shown in Figure 5-23 (1), O), and (k)). Ultimately, the principle crack 
advancement appears to grow along the interfaces located until reaching the top layer. In contrast, 
a second crack pattern tends to branch from broken interfaces of middle layers, leading to shearing 
a few filaments (shown in Figure 5-23 (1), (m), and (n)). The boundary effect on the fracture 
process is likely the primary source of filament cracking and, thus, this secondary crack pattern is 


not considered to be a part of the damage mechanism. 
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Figure 5-23. Mechanical response of Bouligand architecture using the three-point flexural test: 
a) stress-displacement for unrestrained and restrained printed specimens; b) specific modulus of 
rupture (MOR); c) schematic visualization of twisting crack pattern evolution in unrestrained 
and restrained Bouligand architecture with pitch angle y = 3°. 
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5.4.5 Concluding remarks 


In summary, the results presented demonstrate two novel damage mechanisms that can allow 
hcp elements to attain damage-resistant properties when they are subject to restrained shrinkage 
loading. These damage mechanisms are two folds: 1) damage localization mechanism as the first 
layer of protection of hcp elements where the microcracks propagated across the filaments of the 
first layer (1.e., restrained filaments) deflect into the weak interfacial region fabricated between the 
first and second layers rather than deflecting straight into the filaments of the second layer; 2) 
damage delocalization mechanism as the second layer of filaments where an interplay between 
weak interface properties and bio-inspired Bouligand architecture leads to a crack twisting 
mechanism and delocalizes damage spread throughout the specimen. For analysis of the former 
damage mechanism, the length of microcracking advancement along the specimen height is 
measured using three mechanical testing approaches and a non-destructive X-ray uCT imaging 
technique. It is found that the presence of heterogeneous interface in 3D printed hcp elements 
improves the reduction of MOR for restrained specimens by ~ 15% when compared to their cast 
counterparts. The measurement of maximum crack propagation length yields to values of 1 mm 
and 2 mm obtained from the global and local compliance analysis of notched restrained 3D printed 
and cast hcp elements, respectively. The analysis of microstructure characteristics of hcp elements 
reveals that the maximum length of crack propagated along the specimen height is 939 um and 
1,580 um (with a resolution of 15 um) for restrained 3D printed and cast hcp elements, 
complementing the results revealed by mechanical tests. The latter damage mechanism shows a 
further improvement of MOR of Bouligand architectures by ~ 25% when compared to their cast 


counterparts., exhibiting the controlled spread of damage without compromising the strength. 
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6. SUMMARY, CONCLUSION 


The main goal of this study is to understand the underlying physical mechanisms of three 
common durability-related issues of internal frost damage, restrained shrinkage damage, and 
aggregate susceptibility to fracture. We thoroughly invested several topics that are important to 
understand the damage mechanisms, including the role of curvature in governing the freezing of 
pore solution and green phase change materials, the source of competing fracture mechanism in 
heterogeneous mortar and concrete, and anisotropic behavior of 3D printed hcp elements dictated 
by the structural orientation of filaments and the interfacial regions between them on alleviating 
the damage associated with shrinkage cracking. 

First, a thermodynamic-based computational heat transfer model was developed to thoroughly 
investigate the mechanisms that are responsible for macroscopic freeze-thaw behavior of air- 
entrained mortar specimens. The computational results were compared to the experimental ones 
obtained for mortar specimens saturated with pore solution using a low-temperature LGCC. 
Effective medium theory and rule mixture homogenization techniques were utilized to estimate 
the effective thermal properties of composite mortar, where the phase transition of in-pore solution 
at nano and micro levels were considered. The role of curvature, owing to a broad range of pore 
sizes, was considered in calculating the volume fraction of freezable pore solution exposed to 
freezing/thawing cycles using measured absorption—desorption isotherms. Our analysis of the 
synergetic effect between the role of curvature and undercooling before ice nucleation onset 
revealed that the consideration of pore size distribution could reasonably be neglected during the 
freezing process due to undercooling. Conversely, the melting of formed ice indicated a gradual 
process as the temperature increases in both the experimental data and the model with a continuous 
pore size distribution. 

Next, the effect of paraffin oil (sustainable phase change material) incorporation as a TES 
agent on enhancing the macroscopic freeze-thaw behavior of cementitious composites is studied. 
The computational results were compared to the experimental data for mortar specimen tested in 
the LGCC (small-scale) test with a temperature range of —40 °C to +24 °C and concrete slabs 
(large-scale) tested with a temperature range of —10 °C to +10 °C. The numerical model was 
validated by measured data from LGCC experimental data, where we leveraged our understanding 


of the role of curvature combined with the in-pore solution solidification in porous composites. It 
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was concluded that the consideration of a gradual phase transition of PCM with a narrow range of 
transition temperature better predicted the thermal response (temperature profile and associated 
heat flow) of the specimen, while the LWA pore sizes had a negligible effect. Analysis of the slab 
test also showed that the time of freezing at which temperature reduced below 0 °C was decreased 
roughly by 9 h when compared to a slab without PCM, which implied the improvement in the 
freeze-thaw performance of slab. The numerical model was also used to investigate the PCM 
effectiveness on reducing the impact of freeze-thaw cycling within concrete pavements located in 
different regions of the United States. It was observed that the transition temperatures of PCM 
during freezing and melting events are important properties that govern the effectiveness of PCM 
at reducing the impact of freeze-thaw cycles. In other words, they help control how quickly the 
temperature of the pavement drops below the freezing point of pore solution within concrete pores. 
In 122 out of 210 cities investigated, a dose of 170 (kg/m3) of FTO-MT1 PCM would reduce the 
freezing time and depth of concrete pavement by at least 10%, primarily throughout the East, the 
Central, and Pacific Northwest regions. The model results also showed that low transition 
temperatures were more effective, where pavements located in an additional 13 and 44 experienced 
a reduction in freezing time and depth of pavements by at least 10%. These preliminary results 
suggested that the use of PCM-LWA composites as a thermal energy storage system shows 
promise in improving the freeze-thaw performance of concrete pavements and merits further, more 
detailed investigations. 

We, then, studied the key mesostructure morphological and mechanical properties of three- 
phase heterogeneous mortar and concrete to understand better the early-age tensile stiffness and 
strength development, and evolution of competing fracture mechanism in cementitious 
composites. Experimental results indicated that a linear relationship existed between the strength 
and degree of hydration for paste specimens. However, a bilinear relationship is observed as the 
aggregate was added to the paste (1.e., mortar and concrete). The rate of development of tensile 
stiffness contrasts with that of strength, as the fracture of aggregates does not affect the stiffness 
as a function of the degree of hydration. A two-step homogenization model multiscale method 
consisted of continuous and discontinuous homogenization schemes for multiscale modeling of 
failure of quasi-brittle mortar (Level II) and concrete (Level III) having heterogeneous 
mesostructures were used. For Level II mortar, the continuous homogenization scheme has been 


used for quantification of pre-peak and peak stress response of mortar. At the same time, the 
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behavior of a macroscopic cohesive crack is defined upon microscopic localization using the 
discontinuous homogenization scheme. The comparison between numerical solutions and 
experimental observations obtained for tensile stiffness and strength development at discrete early 
ages showed that the variations are less than 5%. With a discontinuous homogenization scheme, 
the fracture energy associated with the amount of energy dissipation over the localized band has 
been estimated to be approximately 15.7 + 3.4 (J/m7). We observed that the pre-peak response of 
mortar is strongly sensitive to the tensile stiffness of their constituents. In contrast, the differences 
between constituents’ strengths lead to govern the micro-crack coalescence and macro-crack 
pattern. We validated the hypothesis of this work by showing that the influence of aggregate 
fracture on the knee point is primarily controlled by the competition between the energy dissipation 
mechanism of the cement paste at very early ages and the energy dissipation mechanism of the 
aggregate phase at later ages. For Level HI concrete, we observed that an increase in area fraction 
of aggregate in two-dimensional samples alleviates the 2D effect, as the variation between 
predicted numerical results and experimental data decreased to less than 5%. Hence, the images 
from the cross-section of the concrete beam tested in the experiment are used to obtain the real 2D 
RVE of concrete using a computed tomography approach. With a continuous homogenization 
scheme, the uncertainty of results obtained from the predictive model for early-age tensile stiffness 
and tensile strength development was less than 10% and 6%, respectively, that indicates an 
outstanding agreement. Similar to mortar, we observed that by understanding the composite nature 
of early age concrete, tensile stiffness and strength development were captured. It was observed 
that the tensile stiffness of mortar and concrete increases at a very similar rate, while the presence 
of coarse aggregates influences the rate of strength development. It was shown that the hydration 
of the paste phase describes the strength gain at an early age; however, it is primarily governed by 
aggregates after a certain degree of hydration is reached. Aggregates do not influence the tensile 
stiffness development similar to how they influence the flexural strength development. 

Finally, we studied the synergetic effect between the weak interface and strong filament 
orientation on enhancing the shrinkage-induced damage tolerant properties of 3D printed 
cementitious systems. As a first step, we elucidated a novel approach detailing the progress in 
creating damage-resistant hardened cement paste (hcp) elements, using additive manufacturing 
(AM). The elements can achieve on-demand strength enhancement when subject to restrained 


shrinkage loading, which is a primary source of the durability-related issue. We incorporate two 
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levels of the structural feature as the two levels of protection into our design. We show that adding 
both first and second levels can greatly enhance mechanical performance (MOR) by 7% compared 
to the design with only the first level of protection and 25% compared to the stiff bulk material. 
Moreover, a detailed FEM analysis for the mechanical response of printed hcp elements with 
lamellar and Bouligand architectures are presented to reveal the damage delocalization mechanism 
evidenced at the layered interfaces. These results and analyses demonstrate proof of concept that 
fabrication of AM-induced weak interfacial regions in conjunction with carefully designed 
structural orientation of bulk filaments offers an unmatched solution for the design of damage- 
tolerant infrastructure elements. Although this study examined the roles of interface and 
architecture in damage resistance of hcp elements under restrained shrinkage loading specifically, 
the methodology presented in this work can be tailored and adapted to gain understanding and 
insight on damage resistance characteristics of other biological materials. Integrating additive 
manufacturing into computational modeling can open a wide new avenue of unprecedented 


application not possible before. 
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